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FOREWORD 


This  is  the  final  report  on  a  study  entitled  "An  Investigation  of  Thick 
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Mechanics,  and  Project  Director. 

Section  II,  an  analysis  of  and  computer  program  for  the  computation  of 
stresses  and  displacements  in  statically  loaded  thick  rings  based  upon  the 
complex  variable  technique  of  the  theory  of  elasticity,  is  the  work  of  Dr. 

Edwin  R.  Chubbuck,  Professor  of  Engineering  Mechanics. 

Dr.  Robert  L.  Thoms,  Associate  Professor  of  Engineering  Mechanics,  pre¬ 
sents  in  Section  III  an  analysis  of  capped  thick  cylindrical  shells  under 
radially  symmetrical  static  loads  using  finite-difference  methods  and  South- 
well  stress  functions. 

Mr.  David  McGill  served  as  checker  and  did  a  major  portion  of  the  pro¬ 
gramming. 
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ABSTRACT 


This  report  describes  three  approaches  to  the  problem  of  predicting  stresses 
and  displacements  in  thick  cylindrical  shells.  Section  I  is  an  analysis  of 
a  ring  or  segment  of  an  infinitely  long  thick  cylindrical  shell  based  upon 
the  simplifying  assumptions  of  the  Winkler  curred  beam  theory.  Dynamic  loading 
of  thick  rings  is  treated  in  Chapter  2  of  Section  I.  Section  II  consists  of 
a  static  analysis  of  the  thick-vailed  circular  cylinder  (or  ring)  by  the  elas¬ 
ticity  approach  developed  by  N.  I.  Muskhelisvili.  Shear  and  radial  stresses 
on  the  inner  boundary,  outer  boundary,  or  both  boundaries  constitute  the 
loading.  A  rather  collate  theoretical  development  is  followed  by  a  computer 
program  and  instructions  for  its  uae.  Section  III  presents  an  analysis  of 
static  stresses  in  axially  loaded  thick-valled  cylinders  with  end  caps.  Ihls 
axlsyametrlc  elasticity  problem  is  solve'  by  finite  difference  techniques  and 
Southwell  stress  functions.  Cylinders  v  -h  one  end  closed  by  either  a  flat 
or  a  hemispherical  cap  are  analysed  and  .  example  worked  for  each  case.  Cylin¬ 
ders  with  both  ends  capped  are  analysed  in  the  final  portion  of  the  report. 
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SECTION  I 


A  STUDY  OF  THE  BEHAVIOR 
OF  THICK  CYLINDRICAL  SHELLS 
USING  WINKLER  CURVED- BEAM  THEORY 

CHAPTER  1 


STATIC  LOADING  OF 
THICK  CIRCULAR  ARCHES  AND  SHELLS 


Note:  The  analysis  which  follows  is  written  to  pertain 

directly  to  a  circular  arch  with  a  plane  of  symme- 
try  and  loaded  symmetrically  such  that  it  deforms 
in  its  original  plane.  The  analysis  is  immediately 
adaptable  to  a  unit  length  of  a  long  cylindrical 
shell  undergoing  plane  strain. 


General 


In  the  analysis  of  curved  beams  or  segments  of  long  cylin¬ 
drical  shells,  simplifying  assumptions  as  to  the  geometry  of 
deformation  are  ordinarily  made.  These  assumptions  allow  one  to 
obtain  relatively  simple  expressions  for  the  stresses  and  de¬ 
formations  resulting  from  loads.  For  design  purposes  the  theory 
has  the  advantage  that  the  significant  quantities,  normal  force, 
shearing  force  and  bending  moment  are  easily  computed. 

The  theory,  commonly  referred  to  as  the  Winkler  curved-beam 
theory,  is  well  known.  It  predicts  tangential  stresses  with 
remarkable  accuracy  except  in  the  neighborhoods  of  concentrated 
loads1  .  There  the  more  powerful  methods  of  the  theory  of  elas¬ 
ticity  are  necessary  to  determine  the  stresses.  The  advantages 
in  the  use  of  the  theory  are  its  simplicity  and  immediate  appli¬ 
cability  to  design;  one  disadvantage  is  that,  as  in  straight- 
beam  theory,  radial  stresses  are  completely  disregarded. 

It  is  the  purpose  of  this  discussion  to  develope  and  present 
in  a  usable  form  the  portions  of  the  theory  of  interest  to  the 
designer. 

We  consider  first  a  segment  of  a  ring  with  constant  curvature 
which  is  symmetrical  with  respect  to  its  center  plane  (x  *  0)  in 
Fig.  1) . 


1 


Figure  1.  Section  of  curved  beam  showing 
coordinates  and  dimensions. 

The  usual  assumptions  of  beam  theory  are  made.  Referring 
to  Fig.  1.  these  are: 

(o)  Hooke's  law  is  valid. 

(b)  ar  is  negligible. 

(c)  The  effect  of  radial  strain  upon  the  strain  due  to 
bending  is  negligible. 

(d)  Sections  normal  to  the  original  centroidal  surface  before 
deformation  are  normal  to  the  new  centroidal  surface 
after  deformation. 

(e)  Displacements  are  small  and  in  the  plane  of  the  arch. 


The  Tangential  Strain 

The  tangential  and  radial  displacement  components  of  particles 
on  the  centroidal  surface  are  denoted  by  v  and  w  respectively,  w 


2 


is  positive  in  the  direction  of  positive  z  (outward)  and  v  is 
positive  in  the  direction  of  increasing  8.  It  can  be  shown8  that 
the  assumptions  lead  to  the  following  expressions  for  v(z)  and 

e  ,  the  tangential  displacement  of  a  particle  and  the  tangential 

0 

strain  of  an  element  at  a  distance  z  from  the  centroidal  surface: 

.  .  v(a+z)  z  dw 

v  (z)  =  — » - u  -  —  — 

a  a  d0 

(1-1) 

_1  dv  z  da  w  w 

9  “  a  d0  a (a+z)  d0a  a+z 

The  radial  displacement,  w(z)  ,  of  a  particle  a  distance  z  from 
the  centroidal  surface  is  taken  equal  to  the  radial  displacement 
of  the  corresponding  particle  on  the  centroidal  surface;  i.e. 
w(z)  *  w. 


The  Strain  Energy  Expression 


The  strain  energy  per  unit  volume  is  *so  e  ,  or,  for  the  arch, 
Ee»  6  9 

-  °  in  which  E  is  Young's  modulus.  (For  the  plane  strain  prob- 

0  e  Eea 

lem,  the  strain  energy  density  is  or  2  j  in  which  u  is 

Poisson's  ratio.  The  strain  energy,  V,  of  a  portion  of  the  arch 
which  subtends  the  angle  a  is 


V 


(a+z)  dAd0 


(1-2) 


Substituting  the  expression  for  e  from  Eq.  (1-1)  into  Eq.  (1-2) 

0 

and  integrating  over  the  area, 


de 


(1-3) 


In  this  expression  A  is  the  area  of  the  beam  cross  section  and 
Z,  the  section  constant  of  curved  beam  theory,  is  defined  by  the 
equivalent  expressions 


_1  f  zdA  =  If  zadA 

A  j.a+z  '  ’  Aa  '  a+z 

A 


Z 


Si  C  dA 
A  J  a+z 


1 


(1-4) 
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The  section  constant  Z  is  dimensionless  and  very  small.  For 
rings  with  geometry  such  that  a  »  z,  the  z  in  the  denominator 
of  the  second  of  equation  (1-4)  can  be  disregarded  in  comparison 
to  a  and 


Z  «  I/aa A, 

in  which  I  is  the  moment  of  inertia  of  the  cross  section  with 
respect  to  the  x  axis.  Evaluation  of  Z  for  a  curved  beam  of  rec¬ 
tangular  cross  section  and  a  mean  radius  to  thickness  ratio  of  3 
reveals  that  Z  %  .009.  As  the  mean  radius  to  thickness  ratio 
increases  Z  approaches  zero. 

Equation  (1-3)  gives  the  strain  energy  of  a  curved  beam  seg¬ 
ment  as  a  function  of  its  configuration  as  defined  by  the  dis¬ 
placement  components  v  and  w  of  particles  on  the  centroidal 
surface.  These  displacements  are  of  course  functions  of  0;  i.e., 
v  =  v(9)  and  w  =  w(9). 

For  complete  rings  and  statically  indeterminate  arches  it 
is  convenient  to  have  the  strain  energy  expressed  as  an  integral 
in  terms  of  the  normal  force  N(0)  and  the  bending  moment  M(0) . 

In  this  form  Castigliano' s  theorem  may  be  used  to  evaluate  re¬ 
dundant  forces  and  moments  and  to  compute  deflections. 

The  normal  force  N(0)is  defined  to  be  positive  when  it  is 
tensile  and  the  bending  moment  M(9)  is  taken  as  positive  when  it 
tends  to  straighten  the  arch.  Thus 

N<8>  -  5  o  dA  =  E  \  e0dA  =  “  +  w  +  z<w  +  fgft  ]  f1’5' 


M  (9 ) 


o  zdA  =  - 

9 


s 


da  w 

eQzdA  =  AEZ  (w  +  rjjpr) 


(1-6) 


Expressing  Eq. (1-3)  in  terms  of  N(9)  and  M(9)  by  means  of  Eqs. (1-5) 
and  (1-6) , 


'  ■  *  -  W  • 


(1-7) 


As  the  radius,  a,  becomes  infinite  this  expression  becomes 


„  l  r~  «■  A  M*  . 

V  *  2  j _  AE  +  El  d*' 
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the  usual  strain  energy  expression  for  a  loaded  straight  beam. 

The  straight  beam  energy  expression  is  commonly  used  in 
arch  analysis;  equation  (1-7)  is  unquestionably  more  accurate. 

Its  use  is  simplified  by  the  fact  that  the  quantity  N(9)  -  ^ 

is  independent  of  9  in  regions  where  there  are  no  distributed 
tangential  loads.  In  the  regions  between  concentrated  tangential 

loads  the  quantity  N(9)  -  is  a  constant.  This  fact  may  be 

verified  by  statics. 


The  Equations  of  Equilibrium 


We  cons  ;.der  now  the  portion  of  a  curved  beam  shown  below  in 

Fig.  2.  It  is  loaded  at  its  outer  surface  by  a  radial  load  q (9 ) , 

a  tangential  load  t(9)  and  by  the  concentrated  forces  and  moment 

N  ,  Q  and  M  at  its  free  end.  The  distributed  loads  are  expressed 
e  e  e 

in  units  of  force  per  unit  of  arc.  The  strain  energy  V  and  the 
potential  energy  of  the  external  loads,  Q,  are  given  below.  The 
total  potential  energy  of  the  system  is  designated  by  U  and  U  = 

V  +  Q.  The  equilibrium  equations  are  obtained  from  the  principle 
of  stationary  potential  energy,  6U  *  0. 
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'-■es* 


+  Z(w  +  |gf)8jd0 
a  d0t}R»d0  "  Ne  V(a)  *  Qe  w(a) 

u  -  v  +  n 


Q  =  -  ^  jqw  +  - 


-M 

e  a  d0 


6 -a 


The  tangential  displacement  v  is  given  a  virtual  displace¬ 
ment  6v  =  €  rj  ( 9 ) ,  consistent  with  the  geometrical  boundary  condi¬ 
tions,  and  the  corresponding  increment  in  U  is  computed. 


U  +  AU 


ffy,ir€T’'+w)s  +  z(w+^)s}d0- 

f  fr^  dS  -  Ne (v+e  ">[  -Qe  w  (a)  -Me  7  ae  I 

0*a  'e*a 


i 


AU  *  6U  +  0  (ca ) 

6u  -  *r S0(d9  +w)T,'de  -  «So  -r-tT^9  -  vn<a) 

In  this  expression  and  what  follows  differentiation  of  ri  with 
respect  to  0  is  denoted  by  tj'. 


In  order  for  6U  to  vanish 

+  iK)  +  At  .  o 

a  'dea  d0  a 

(1-8) 

„  AE,dv  A  A 

e  a  'd0 

(1-9) 

The  radial  displacement  w  is  now  varied  by  6w  ■  c r|  (0 ) 


-  6  - 


and 


U  +  AU  *  fa  5  {*d^  +  W+  ex])9+  Z(w+  €T1  + d0^  +  GT1*)a} 

.a 


d2  w 


de 


-  $  {q(w+er,)  +  h*  -  f*  (ds  +  «„')}* 


dQ 


-  N  v(a)  -  Q  (w+eri) 

C 


|9  ^OL 


Me,dw  , . 

-  “ (d5"  +  €T1  > 


e=a 


6U  = 


AEe 


L{(*T  +  w)ri  +  2<w  +  +  z<«  +  |jr>n'}de 


ra  r  1 

-  ^  {q€Ti  -  ^-eTi'j-Rode  “  Qecri(a)  -  ‘jneV  (a) 


For  6U  to  vanish 


d4  w  „da  w 


+  2—  +  w  +  +  W>-  gR"a  cR"  dt  _  0 

de4  ‘de3  +  z  de  +  aez  aez  de  0 


M  = 


Q  = 


da  w 

AEZ  (w  ♦  fjft 


8=a 


ctRo 

a 


0  =a 


AEZ  .dw  d3  w 
a  dp  +  dp* 


9=a 


(1-10) 

(1-11) 

(1-12) 


Equations  (1-8)  through  (1-12)  are  thus  the  equilibrium  equa¬ 
tions  and  the  boundary  conditions  for  the  loaded  arch. 

Equation  (1-8)  states  that 


M  ±-,*L  +  w)  =  .  &£ 
a  de  lde  '  a  • 


dv 


If  t  i  0  in  a  region  then  —  +  w  is  a  constant  in  that  region 
From  equations  (1-5)  and  (1-6) 

N,e>  -&&-  -a  (&  +w). 
a  a  ae 

M  (0 ) 

Thus  N(0)  -  — is  a  constant  in  regions  where  there  is  no 
tangential  load. 
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Equation  (1-9)  simply  states  that  the  normal  force  at  the 

end,  N  ,  equals  AE  times  the  strain  of  the  centroidal  surface, 
e 

The  solution  to  any  equilibrium  problem  of  a  ring  or  arch 
so  loaded  thus  may  be  obtained  by  solving  for  the  functions  v 
and  v  from  equations  (1-8)  and  (1-10)  subject  to  the  geometrical 
boundary  conditions  and  the  natural  boundary  conditions  given  by 
equations  (1-9) ,  (1-11) ,  and  (1-12) . 


The  Tangential  Stress 


Once  v  and  w  are  known,  a.  may  be  computed  from  a  =  Ee 

6  0  0 


with  e  given  by  equation  (1-1) 
9 


It  is  convenient,  however,  to 


have  an  expression  for  in  terms  of  the  normal  force  and  bend¬ 
ing  moment  at  any  section,  a.  =  Ee  or 

9  9 


a 


0 


Er . 

z 

d3  w 

J 

— 1 

Lde  a 

a  (a+z) 

de3 

a+zj 

prl  dv  _ 

z 

d3  w 

w 

zw 

ELa  de 

a  (a+z) 

dF 

a+z 

a (a+z) 

ZW  ~[ 

a (a+z) J 


rl  dv  w 

El_7  as  +  7 


Z 

a (a+z) 


(w  + 


From  equations  (1-5)  and  (1-6) 


M  (8 )  z  I 
aZ (a+z) J 


(1-13) 


The  Complete  Ring  with  Radial  Loads 

We  consider  now  a  complete  ring  loaded  at  its  outer  surface 
by  radial  and  tangential  loads  whose  resultant  is  zero.  Castig- 
liano's  theorem  and  the  principle  of  superposition  permit  one  to 
determine  the  normal  and  shear  forces  and  the  bending  moment  at 
any  section  and  to  analyze  deflections. 

We  consider  first  the  ring  shown  below  loaded  by  only  a  con¬ 
centrated  force  P  at  the  angle  a.  P  is  positive  in  what  follows 
when  as  shown. 
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Figure  3.  Complete  ring  fixed  at  0  =  0  and 
with  radial  load  at  6  =  a 

The  ring  is,  of  course,  not  in  equilibrium  under  the  action 
of  P  and,  for  convenience,  we  imagine  a  cross  section  at  0  =  0 
to  be  fixed  by  external  supports.  The  first  problem  considered 
is  that  of  determining  the  redundant  forces  and  moment  N,  Q  and 
M  at  0  =  +0  due  to  P  at  a.  Having  these  quantities  due  to  a 
concentrated  load  at  a,  they  may  be  determined  for  any  system  of 
concentrated  radial  loads  by  superposition.  A  distributed  radial 
load  may  be  treated  by  replacing  it  with  an  equivalent  system  of 
concentrated  loads.  If  the  loading  system  is  self  equilibrating 
the  external  supports  required  at  0  =  0  to  maintain  equilibrium 
vanish. 

The  expressions  for  N (0 )  and  M(0)  at  any  angle  0  due  to  P 
at  a  are  given  below. 

N  (0 )  =  N  Cos  0  -  Q  Sin  0  for  0  <  0  <  a 

M (0 )  =  M  -  Na (1-Cos  0)  -  Qa  Sin  0  for  0  <  0  <;  a 

N (0 )  =  N  Cos  0  -  Q  Sin  0  -  P  Sin(0-a)  for  a  <  0  £  2rr 

M (0 )  =  M  -  Na  (1-Cos  0)  -  Qa  Sin  0  -  Pa  Sin(e-a) 

for  a  *  0  s  2tt 
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The  strain  energy  from  equation  (1-7)  is  thus 


Since  there  is  no  displacement  or  rotation  of  the  cross  section 
upon  which  M,  N  and  Q  act,  it  follows  that 


av  _  av  av 
aM  dN  ao 


Setting  these  derivatives  of  V  equal  to  zero  and  simplifying,  one 
obtains  the  three  influence  functions  for  M,  N  and  Q  given  below. 


M 

Pa 


1) Sin  a  + 


(1-Cos  a) 
2tt(1+Z) 


N 

P 


l)Sin  a 


Q  Sir?  a  ,  a  Sin  2a.  Cos  a 

P  — s;  (n  - 1 +  4  « 


(1-14) 


It  is  to  be  noted  that  N  and  Q  are  completely  independent  of 
the  thickness  of  the  ring.  Since  Z  is  very  small,  M  also  is 
essentially  unaffected  by  the  thickness  to  radius  ratio.  This 
fact  is  important  since  many  solutions  based  on  ordinary  arch 
theory  exist.  They  can  be  used  with  great  accuracy  for  determin¬ 
ing  the  redundant  reactions  due  to  radial  loads. 


Knowing  M,  N  and  Q  one  may  draw  complete  shear  and  bending 
moment  diagrams  for  the  ring.  Equation  (1-13)  will  give  o  at 

n 

any  cross  section.  Elementary  theory  is  ordinarily  used  for  comput 
ing  shearing  stresses. 


To  facilitate  computation,  the  table  of  influence  coefficients 
is  presented  in  Table  1  which  follows. 
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TABLE  I* 


M,  N,  and  Q  at  0  =  +0  Due  to 
P  at  a 


a0 

Q/P 

N/P 

M/Pa 

0 

-1.00000 

.000000 

.000000 

10 

-.985089 

-.168825 

-.166407 

20 

-.941922 

-.323019 

-.313421 

30 

-.873434 

-.458333 

-.437011 

40 

-.783231 

-.571367 

-.534132 

50 

-.675431 

-.659649 

-.602797 

60 

-.554499 

-.721688 

-.642110 

70 

-.425073 

-.756975 

-.652254 

80 

-.291797 

-.765962 

-.634444 

90 

-.159155 

-.750000 

-.590845 

100 

-.031324 

-.711250 

-.524458 

110 

.087957 

-.652564 

-.438975 

120 

.195501 

-.577350 

-.338618 

130 

.288750 

-.489417 

-.227960 

140 

.365836 

-.392815 

-.111740 

150 

.425604 

-.291667 

.005320 

160 

.467617 

-.190011 

.118700 

170 

.492123 

-.091648 

.224244 

180 

.500000 

.000000 

.318310 

190 

.492685 

.082001 

.397892 

200 

.472075 

.152009 

.460721 

210 

.440421 

.208333 

.505321 

220 

.400209 

.249973 

.531048 

230 

.354038 

.276627 

.538085 

240 

.304499 

.288675 

.527408 

250 

.254063 

.287128 

.500717 

260 

.204973 

.273558 

.460350 

270 

.159155 

.250000 

.409155 

280 

.118148 

.218846 

.350364 

290 

.083053 

.182718 

.287439 

300 

.054499 

.144338 

.223915 

310 

.032644 

.106395 

.163247 

320 

.017187 

.071421 

.108656 

330 

.007409 

.041667 

.062989 

340 

.002229 

.019001 

.028599 

350 

.000281 

.004824 

.007241 

360 

.000000 

.000000 

.000000 

♦This  table  was  computed  from  equations  (1-14)  with  Z  set 
equal  to  zero. 
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In  the  computation  of  this  table,  Z  was  set  equal  to  zero 
since  it  is  always  small  in  comparison  to  unity.  Unfortunately, 
when  tangential  loads  are  considered,  the  thickness  to  radius  ratio 
is  of  significance  and  influence  tables  are  necessary  for  each 
thickness  to  radius  ratio. 


The  Complete  Ring  with  Tangential  Loads 

Referring  to  Figure  (3)  and  replacing  the  radial  load  P  at 
a  with  a  concentrated  tangential  load  T  (positive  in  the  direction 
of  increasing  0)  at  a  then 

N  (0 )  =  N  Cos  0  -  Q  Sin  0  for  0  <  0  <  a 

M(e)  =  M-Na (1-Cos  0)  -  Qa  Sin  0  for  0  <  0^a 

N  (9 )  =  N  Cos  0  -  Q  Sin  0  -  T  Cos(0-a)  for  a  <  0  ^  2tt 

M(0)  =  M-Na (1-Cos  9)  -  Qa  Sin  0  +  Tc  +  Ta[ 1-Cos  (0-a) ] 

for  a  ^  9  £  2tt 


The  strain  energy  by  equation  (1-7)  is  thus 

a 


v*.  •> 

2EAV  =  \  |a(N-^)2  +  *“■  [M-Na (1-Cos  0)-Qa  Sin  0]a]>d0 

+  ^  ^a[^N-^-T(l+c/a)  J  +  ~  ^M-Na  (1-Cos  0)-Qa  Sin  9 

a 

(e-a)]a}de 


+  Tc  +  Ta  -  Ta  Cos 


Again,  since  there  is  no  displacement  or  rotation  of  the 
cross  section  upon  which  M,  N  and  Q  act 

av  _  av  _  av 

dM  0N  dQ 

Setting  these  derivatives  of  V  equal  to  zero  and  simplif  ng 


M  1  fZ  Sin  a  c 

57  =  n  L2  (1+Z)  +  “  Sin  a  +  n  Cos  a 

+  f(l-Cos  a)  -nU+f)  +  §J] 

T  =  n  Hi  +  7)Sin  a  +  COB  a] 


(1-15) 
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f  !l-Coa  a)  -  (l-^-)Sin  a 

The  influence  tables  presented  below  were  computed  using 

R-.  D  h 

c  =  y  L  =  — .  Tables  are  presented  for  thickness,  h,  to  mean 
22  111  1 

radius,  a,  ratios  of  *r,  — ,  — ,  and  — .  Linear  interpolation  may 

3  4  5  o 

be  used  with  accuracy  for  intermediate  values. 
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TABLE  II 


M,  N,  and  Q  at  8  =  +0  Due  to 
T  at  a.  h/a  *  1/3 


a0 

Q/T 

N/T 

M/Ta 

0 

.000000 

1.000000 

-.166667 

10 

-.174466 

.994301 

-.167337 

20 

-.345415 

.960066 

-.195712 

30 

-.508086 

.899960 

-.248320 

40 

-.658249 

.817332 

-.321053 

50 

-.792304 

.716071 

-.409341 

60 

-.907369 

.600443 

-.508326 

70 

-1.001323 

.474925 

-.613051 

80 

-1.072837 

.344042 

-.718640 

90 

-1.121362 

.212207 

-.820464 

100 

-1.147098 

.083570 

-.914297 

110 

-1.150939 

-.038105 

-.996452 

120 

-1.134393 

-.149557 

-1.063881 

130 

-1.099485 

-.248110 

-1.114263 

140 

-1.048656 

-.331735 

-1.146046 

150 

-.984637 

-.399078 

-1.158469 

160 

-.910338 

-.449473 

-1.151547 

170 

-.828729 

-.482910 

-1.126030 

180 

-.742723 

-.500000 

-1.083333 

190 

-.655081 

-.501897 

-1.025444 

200 

-.568318 

-.490220 

-.954812 

210 

-.484637 

-.466947 

-.874223 

220 

-.405868 

-.434310 

-.786665 

230 

-.333441 

-.394677 

-.695191 

240 

-.268367 

-.350443 

-.602785 

250 

-.211246 

-.303915 

-.512235 

260 

-.162290 

-.257218 

-.426017 

270 

-.121361 

-.212207 

-.346203 

280 

-.088029 

-.170394 

-.274379 

290 

-.061630 

-.132905 

-.211595 

300 

-.041343 

-.100443 

-.158341 

310 

-.026260 

-.073284 

-.114538 

320 

-.015461 

-.051288 

-.079569 

330 

-.008086 

-.033934 

-.052322 

340 

-.003395 

-.020374 

-.031262 

350 

-.000818 

-.009493 

-.014522 

360 

.000000 

.000000 

.000000 
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TABLE  III 


M,  N,  and  Q  at  0  =  +0  Due  to 
T  at  a.  h/a  =  1/4 


a0 

Q/T 

N/T 

M/Ta 

0 

.000000 

1.000000 

-.125000 

10 

-.174265 

.991998 

-.129244 

20 

-.344615 

.955530 

-.161119 

30 

-.506310 

.893328 

-.217083 

40 

-.655146 

.808807 

-.292961 

50 

-.787567 

.705911 

-.384121 

60 

-.900737 

.588957 

-.485654 

70 

-.992596 

.462462 

-.592562 

80 

-1.061877 

.330981 

-.699936 

90 

-1.108099 

.198944 

-.803129 

100 

-1.131532 

.070509 

-.897909 

110 

-1.133140 

-.050568 

-.980593 

120 

-1.114498 

-.161043 

-1.048154 

130 

-1.077697 

-.258270 

-1.098302 

140 

-1.025233 

-.340260 

-1.129528 

150 

-.959888 

-.405710 

-1.141121 

160 

-.884612 

-.454009 

-1.133158 

170 

-.802405 

-.485213 

-1.106456 

180 

-.716197 

-.500000 

-1.062500 

190 

-.628756 

-.499594 

-1.003352 

200 

-.542592 

-.485684 

-.931534 

210 

-.459888 

-.460316 

-.849904 

220 

-.382445 

-.425785 

-.761517 

230 

-.311653 

-.384517 

-.669485 

240 

-.248473 

-. j Jb957 

-.3/6846 

250 

-.193447 

-.291452 

-.486427 

260 

-.146724 

-.244157 

-.400740 

270 

-.108099 

-.198944 

-.321871 

280 

-.077069 

-.157333 

-.251415 

290 

-.052904 

-.120442 

-.190417 

300 

-.034712 

-.088957 

-.139346 

310 

-.021522 

-.063124 

-.098091 

320 

-.012358 

-.042762 

-.065995 

330 

-.006309 

-.027303 

-.041892 

340 

-.002595 

-.015838 

-.024188 

350 

-.000617 

-.007190 

-.010948 

360 

.000000 

.000000 

.000000 
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TABLE  IV 


M,  N,  and 

Q  at  9  =  +0  Due 

to 

T  at 

a.  h/a  *  1/5 

a0 

Q/T 

N/T 

M/Ta 

0 

.000000 

1.000000 

-.100000 

10 

-.174144 

.990616 

-.106373 

20 

-.344135 

.952809 

-.140333 

30 

-.505243 

.889350 

-.198295 

40 

-.653284 

.803692 

-.276047 

50 

-.784724 

.699815 

-.368919 

60 

-.896758 

.582065 

-.471973 

70 

-.987360 

.454984 

-.580184 

80 

-1.055301 

.323144 

-.688625 

90 

-1.100141 

.190986 

-.792637 

100 

-1.122192 

.062672 

-.887986 

110 

-1.122460 

-.058046 

-.970992 

120 

-1.102562 

-.167935 

-1.038640 

130 

-1.064624 

-.264366 

-1.088656 

140 

-1.011179 

-.345375 

-1.119558 

150 

-.945036 

-.409689 

-1.130667 

160 

-.869177 

-.456730 

-1.122094 

170 

-.786610 

-.486595 

-1.094695 

180 

-.70028? 

-.500000 

-1.050000 

190 

-.612962 

-.498213 

-.990112 

200 

-.527157 

-.482962 

-.917599 

210 

-.445038 

-.456337 

-.835359 

220 

-.368391 

-.420670 

-.746486 

230 

-.298580 

-.378422 

-.654131 

240 

-.236536 

-.332065 

-.561360 

250 

-.182768 

-.283974 

-.471028 

260 

-.137384 

-.236320 

-.385662 

270 

-.100141 

-.190986 

-.307363 

280 

-.070493 

-.149496 

-.237727 

290 

-.047668 

-.112964 

-.177796 

300 

-.030733 

-.082065 

-.128027 

310 

-.018680 

-.057028 

-.088293 

320 

-.010496 

-.037647 

-.057909 

330 

-.005243 

-.023324 

-.035679 

340 

-.002115 

-.013116 

-.019975 

350 

-.000496 

-.005809 

-.008820 

360 

.000000 

.000000 

.000000 
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TABLE  V 


M,  N,  and 

Q  at  8  *  +0  Due 

to 

T  at 

a.  h/a  =  1/8 

a0 

Q/T 

N/T 

M/Ta 

0 

.000000 

1.000000 

-.062500 

10 

-.173963 

.988544 

-.072044 

20 

-.343415 

.948726 

-.109110 

30 

-.503644 

.883381 

-.170051 

40 

-.650492 

.796019 

-.250595 

50 

-.780460 

.690671 

-.346020 

60 

-.890790 

.571728 

-.451341 

70 

-.979506 

.443768 

-.561497 

80 

-1.045437 

.311389 

-.671533 

90 

-1.088204 

.179049 

-.776773 

100 

-1.108183 

.05C917 

-.872977 

110 

-1.106441 

-.069263 

-.956472 

120 

-1.084657 

-.178272 

-1.024258 

130 

-1.045015 

-.273510 

-1.074090 

140 

-.990098 

-.353048 

-1.104523 

150 

-.922764 

-.415657 

-1.114922 

160 

-.846023 

-.460813 

-1.105454 

170 

-.762918 

-.488668 

-1.077033 

180 

-.676408 

-.500000 

-1.031250 

190 

-.589270 

-.496140 

-.970275 

200 

-.504003 

-.478880 

-.896739 

210 

-.422764 

-.450369 

-.813603 

220 

-.347311 

-.412997 

-.724022 

230 

-.278970 

-.369278 

-.631197 

240 

-.218631 

-.321728 

-.538242 

250 

-.166749 

-.272757 

-.448048 

260 

-.123375 

-.224565 

-.363171 

270 

-.088204 

-.179049 

-.285727 

280 

-.060629 

-.137741 

-.217319 

290 

-.039814 

-.101747 

-.158983 

300 

-.024765 

-.071728 

-.111159 

310 

-.014416 

-.047884 

-.073692 

320 

-.007704 

-.029975 

-.045861 

330 

-.003644 

-.017356 

-.026424 

340 

-.001395 

-.009033 

-.013698 

350 

-.000315 

-.003736 

-.005649 

360 

.000000 

.000000 

.000000 
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Ring  with  concentrated  loads 

Referring  to  Table  I,  M,  N  and  Q  due  to  the  load  P  at  90°  are 

M  =  -.5908  Pa 
N  =  -.7500  P 
Q  =  -.1592  P 

The  values  due  to  2T  at  180°  are 

M  *  2  ( .3183)  Ta 
N  =  2 (.0000)  T 
Q  *  2  (.5000)  T 

The  values  due  to  P  at  270°  are 

M  =  +.4092  Pa 
N  =  +.2500  P 
Q  =  +.1592  P 
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Due  to  T  at  90°  from  Table  II 

M  =  -.8205  Ta 
N  =  +.2122  T 
Q  =  -1.1214  T 

Due  to  -T  at  270°  from  Table  II, 

M  =  -.3462  (-Ta) 

N  =  -.2122  (-T) 

Q  =  -.1214  (— T) 

By  superposition 

M  =  -.1816  Pa  +.1623  Ta 
N  =  - . 5P  +  .4144  T 
Q  =  0. 

B.  As  a  second  example,  we  consider  the  ring  loaded  as  shown 
below. 


Ring  with  distributed  radial  load. 
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This  distributed  radial  load  could  be  replaced  by  an  equivalent 
set  of  concentrated  loads  and  Table  I  could  be  used.  In  this  case, 
however,  the  loading  function  is  simple  and  M,  N  and  Q  may  be 
easily  computed  from  equations  (1-14)  by  integration.  For  example, 
from  equation  (1-14)  ,  neglecting  Z, 

M  - pa  - 1)sin  a  +  l=*r ~  ] 


Replacing  P  by  W  Sin  a  da  from  0  to  n  and  by  -W  Sin  c.  Ro  da 
from  n  to  2rr 


dM  =  W  Sina 


l)Sin  a 


+ 


1-Cos  a~ 
2n 


Thus 


M 


=  WRo  a  ^  |(^“  -  l)Sin  a  +  1  -■C"  ^-a  jsin  a  da 

2  TT  .  _ 

-  WR^.a  ^  [(~  -  1)  Sin  a  +  ~  ^--jsin  a  da 

a 


Integrating 

M  =  -WRq  a  (rt/4  -  ~)  . 

N  and  Q  may  be  computed  in  a  similar  fashion. 

Radial  Displacements  due  to  Radial  Loads 

We  consider  the  complete  ring  shown  below.  Fig.  4,  with  a 
radial  load  P,  taken  positive  inward  as  shown,  at  a,  and  a  concen¬ 
trated  radial  load  R  at  3 . 


Figure  4.  Ring  with  radial  loads  at  8=a  and  e=6. 
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M,  N  and  Q  are,  as  before,  the  redundant  reactions  at  0  =  +0. 
Then  for  8  £  a 

N(9)  =  N  Cos  0  -  Q  Sin  0  for  0  <  0  <  8 

M (0 )  =  M  -  Na (1-Cos  0)  -  Qa  Sin  0  for  0  <  0  <  0 

N(9)  =  N  Cos  0  -  Q  Sin  0  -  R  Sin(0-8)  for  8  <  9  <  a 

M (0 )  =  M  -  Na (1-Cos  0)  -  Qa  Sin  0  -  Ra  Sin(0-8)  for 

8  £  9  *  a 

N (0 )  =  N  Cos  9  -  Q  Sin  0  -  R  Sin (0-8)  -  P  Sin(9-a)  for 

a  <  9  s.  2tt 

M (9 )  =  M  -  Na(l-Cos  9)  -  Qa  Sin  0  -  Ra  Sin(9-8)  -  Pa  Sin(q-a) 
for  a  <  0  <;  2rr 


The  strain  energy  is  then  given  by 

2EAV  =  \  ja(N-^)  +  [M-Na  (1-Cos  9)-Qa  Sin  9]  ]• 


d0 


+  ^  |a(N-^)a4^-  M-Na  (1-Cos  9)  -Qa  Sin  9-Ra  Sin(0-8)]  }d9 


8 
2tt 


+  ^  {a(N-^)a+  J^M-Na (1-Cos  0)-Qa  Sin  Q-Ra  Sin(0-8) 

a 


-  Pa  Sin(9-a)J  J-d0 


dV 

The  radial  deflection,  w,  at  8  due  to  P  at  a  equals  —— ) 

.  dV.  8  R=0 

dR  R=0 

Performing  this  differentiation  and  integration  and  solving 
dV 

for  rrr  with  R  set  equal  to  zero, 

c>R 

EAZ  w(B,a)  =  (M-Na) (Cos  0-1)  +  (n-0/2)Na  Sin  0 

+  (tt— 8/2)  Qa  Cos  8  +  ^  Sin  B 

Pa 

+  (n-a/2)  Pa  Cos(0-a)  +  —  Cos  8  Sin  a. 
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Substituting  for  M,  N  and  Q  their  values  from  equations  (1-14) 
and  simplifying 


EAZ 

Pa 


w(0  ,a) 


(C0»  g-1) (l-Co«  B)  fi_,  a. 

2tr  (1+Z)  2tt  '  2' 


Cos (0-a) 


+ 


1_ 

4n 


Cos  0  Sin  a 


Sin  B  Sin  a 
4tt 


Sin  8  Cos  a 
2tt 


(1-16) 

Equation  (1-16)  is  valid  for  the  radial  deflection  at  0  due 
to  a  radial  load,  P,  at  a  with  0  £  a.  The  computations  were  re¬ 
peated  for  the  case  in  which  0  ^  a  and  in  this  case 


EAZ 

Pa 


w(0,a) 


(Cos  0-1)  (1-Cos  g) 
2tt  (1+Z) 


+  ^■(tt-0/2)  Cos  (0-a) 


+  — —  Cos 
4n 


Sin  0  Sin  a  , 

Sin  0  -  - -  -  (rr-0/2) 


Sin  a  Cos  a 


2tt 


(1-17) 

Equations  (1-16)  and  (1-17)  are  the  complete  influence  function 
for  the  radial  displacement  due  to  radial  loads  for  the  entire  ring. 
It  is  to  be  noted  that  w(0,a)  =  w(a,0);  i.e.,  if  a  and  0  are  inter¬ 
changed  in  either  equation  (1-16)  or  (1-17)  the  equations  become 
identical.  Also,  from  Fig.  4  it  is  obvious  that  w(0,a)  = 
w(2tt-0 , 2rr-a)  . 

From  equation  (1-17)  a  complete  set  of  influence  coefficients 
with  10° increments  and  with  Z  on  the  right  hand  side  set  equal  to 
zero  has  been  computed.  By  omitting  0°  and  360° ,  this  set  is  a 
35  x  35  array,  symmetrical  with  respect  to  both  diagonals. 

The  following  program,  written  in  Fortran  II,  will  generate 
this  matrix.  The  deflection  due  to  any  radial  loading  may  be  com¬ 
puted  by  constructing  a  column  loading  matrix  and  multiplying  it  by 
the  influence  matrix. 


22 


DIMENSION  P  (35, 35) , JJ (35) 

DO  1  1=10, 180, 10 

K=I/10 

ALPHA=I 

ALPHA= ALPHA* .017453293 
SA=SINF (ALPHA) 

CA=COSF (ALPHA) 

IJ=360-I 

DO  1  J=I,IJ,10 

L=J/10 

BETA=J 

BETA=BETA*. 017453293 
SB=SINF (BETA) 

C3=C0SF (BETA) 

COSAB=COSF  (ALPHA-BETA) 

1  OP (K, L) =-  (CB-1 . )  * (CA-1 . ) * . 15915494+ ( . 5-BETA* . 079577471)  * (ALPHA* 
1C0SAB-SA*CB) + (ALPHA*CA-SA) *SB* .079577471 

DO  2  1=1,17 
IJ=35-I 
K=36-I 
DO  2  J=I,IJ 
M=36-J 

2  P  (M,K)  =P  (I ,  J) 

DO  3  1=2,35 
IJ=I-1 

DO  3  J=1,IJ 

3  P(I#J)=P(J#I) 

4  FORMAT (13, 1HO, 6X, 7F10. 6) 

5  FORMAT  (5H . ) 

6  FORMAT  (15HBETA  .  ALPHA  ,  6  (13 , 7X)  ,  13) 

7  FORMAT  (12H  ) 

Ml=l 

M2=7 

DO  9  N=1 , 5 
PUNCH  5 
DO  8  I=M1,M2 

8  JJ(I)=I*10 

PUNCH  6,  (JJ(I) ,I=M1,M2) 

PUNCH  7 

PUNCH  4,  (J,  (P(I,J) ,I=M1,M2)  ,J=1,35) 

Ml=M2+l 

9  M2=M2+7 
END 
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Radial  Displacements  due  to  Tangential  Loads 

We  now  consider  the  complete  ring  loaded  as  shown  below  in 
Fig.  5.  It  carries  a  concentrated  tangential  load,  T,  at  the 
angle  a  and  a  concentrated  radial  load,  R,  at  the  angle  0  with 
a  ^  8.  The  radial  deflection  at  8  due  to  T  at  a  equals  dV/dR 
with  R  set  equal  to  0. 


T 


Figure  5.  Ring  with  radial  load  at  8  and 
tangential  load  at  a. 

The  values  of  M(9)  and  N(0)  at  any  angle  0  are  as  given  below. 

N(0)  =  N  Cos  0  -  Q  Sin  0  for  0  <  0  <  8 

M (0 )  ■  M  -  Na (1-Cos  9)  -  Qa  Sin  0  fou  0  <  0  £  0 

N(0)  *  N  Cos  0  -  Q  Sin  0  -  R  Sin(0-0)  for  8  <  0  <  a 

M (0 )  *  M  -  Na(l-Cos  0)  -  Qa  Sin  9  -  Ra  Sin(9-0)  for 

8  £  9  *  a 
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N  (9 )  =  N  Cos  0  -  Q  Sin  9  -  R  Sin  (9-8)  -  T  Cos(0-a)  for 
a  <  0  *  2tt 

M (0 )  =  M  -  Na(l-Cos0)  -  Qa  Sin  9  -  Ra  Sin(0-8)  +  Tc 
+  Ta[  1-Cos  (9-a)  ]  for  a  *  0  £  2tt. 

Then 

2EAV  =  \  |a(N-^)2  +  ^M-Na (1-Cos  0)  -  Qa  Sin  ©J  j-d0 

+  \  {a(M--^)a+  ^  [M-Na (1-Cos  0 )  -  Qa  Sin  0  -  Ra  Sin  (0-8)  J* }d0 
a 

C  f  r  M  Cl2  1  r 

+  \  |al  N-  —  -  T  (1+  -)\  +  —  I  M-Na  (1-Cos  0)  -  Qa  Sin  0 
a 

-  Ra  Sin  (0-8)  +  Tc  +  Ta  -  Ta  Cos(0-a)J  j>d0 

Differentiating  the  above  expression  for  V  with  respect  to 
R,  setting  R  =  0,  substituting  the  values  of  M,  N  and  Q  from 
equations  (1-15) ,  and  simplifying,  there  results. 


fr  1$  -  fr  v(8'a)  =  <1+f> <1-“/2t’>  f1-008  b> 

Sin  a  rSin  8  Sin  a,l 

+  2n7I+zT  (1-Cos  +  „ - <2  +  c/») 

+  (l-a/2n) Sin (8-a) -  (1+c/a) (1-Cos  a)C°°  Pj 
-  (1+c/a)  ^  (1-Cos  a)+Cos(8-a)  -  Cos  b] 


(1-18) 


-  (tt-ol/2)  Sin  (8-a) 


+  7—  Sin  8  Sin  a. 
4rr 


Equation  (1-18)  is  valid  for  computation  of  radial  deflections 
at  the  angle  8  resulting  from  a  tangential  load,  T.  at  the  angle  a 
with  a  *  8 .  The  problem  was  reworked  with  the  loading  given  in 
Fig.  5  but  with  0  ^  a.  The  final  result  only  is  presented  here. 
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EAZ 

Ta 


w(0,a)  »  (Cos  e-l)[2n(1+ a)  -  2tt  (1+Z)  ] 

-  (n  -|)[C-°-°  ^  (1+c/a)  (1-Cos  a)  -  —  Sin(a-S) 

Sin  a  Sin  8  ,1  .  ,  „  i  Sin  8  _ 

-  - - - ^(j  +c/a)  I-  -  (1-t-c/a)  (1-Cos  a) 

+  7—  Sin  8  Sin  a 

4n 


(1-19) 


Equations  (1-18)  and  (1-19)  represent  the  entire  influence 
function  for  a  complete  ring  carrying  tangential  loads  on  its 
outer  surface.  In  this  case  w((3,a)  /  w(a,9)  since  the  radial 
displacement  at  8  due  to  a  tangential  load  at  a  is  not  equal  to 
the  radial  displacement  at  a  due  to  a  tangential  load  at  8* 

There  is,  however,  a  different  type  of  symmetry.  Consideration 
of  Fig.  5  reveals  that  w(0,a)  =  -w(2n-8*  2n-a)  . 

The  following  program,  written  in  Fortran  II,  will  generate 
a  complete  set  of  influence  coefficients  for  radial  displacements 
due  to  tangential  loads  for  10°  increments.  In  the  program,  c  was 
set  equal  to  h/2  and  h/a  may  be  assigned  any  value.  The  matrix 
is  antisymmetrical  with  respect  to  its  center  element,  a  fact 
which  was  used  to  shorten  the  program.  The  value  of  the  radial 
deflection  at  each  10°  interval  may  be  obtained  by  constructing 
a  loading  column  matrix  and  multiplying  it  by  the  influence 
matrix. 
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DIMENSION  P (35, 35) , JJ(35) 

1  READ  8,  TOR 
PUNCH  9, TOR 
P (18, 18) =0 . 

CONST=l . + . 5*TOR 
HALF=(l.+TOR) *. 5 

Z=l./( (l./TOR)*LOGF((2.+TOR)/(2.-TOR) )) 

M=1 

DO  3  1=10,340,10 

IP1=I+10 

GO  TO  4 

2  DO  3  J=IPl , 350 , 10 
GO  TO  5 

3  CONTINUE 
GO  TO  6 

4  K=I/10 
Il=36-K 
ALPHA=I 

ALPHA=ALPHA* . 0 1 74 5 3 2 9 3 
SA=SINF  (ALPHA) 

CA=1 . -COSF (ALPHA) 

ALPH2=ALPHA* .  5 
GO  TO  (2,5) ,M 

5  L=J/10 
I2=36-L 
BE  TA= J 

BETA=BETA* . 017453293 
SB=SINF  (BETA) 

CB=C0SF  (BETA) 

CBl=l . -CB 

CBCA=CB*CA*CONST 

SBSA=SB*SA 

SBMA=SINF  (BETA- ALPHA) 

PMB=3 . 14 15927-BETA* . 5 

OP  (K,L)=  ( (SA*.  5*Z-ALPH2*CONST)  *CBl+  (HALF*SBSA-ALPH2*SBMA-CBCA)  * 
1PMB- (CONST*CA-ALPH2*SA) *SB* . 5) *. 31830989 
P(I1,I2)=-P(K,L) 

GO  TO  (3, 7) ,M 

6  M=2 

DO  7  1=10,170,10 
J=I 

GO  TO  4 

7  CONTINUE 

8  FORMAT (E15 . 8) 
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9  FORMAT  ( 3 7HRATI0  OF  THICKNESS  TO  MEAN  RADIUS  IS  ,F12.8) 

10  FORMAT  (5H . ) 

11  FORMAT  (15HBETA  .  ALPHA  ,  6  (13 , 7X)  ,  13) 

12  FORMAT  (12H  . ) 

13  FORMAT  (13 , 1H0,6X,  7F10.6) 

Ml=l 

M2=7 

DO  15  N-1,5 
PUNCH  10 
DO  14  I=M1 ,M2 

14  JJ(I)=I*10 

PUNCH  11,  (JJ(I) , I=Ml,M2) 

PUNCH  12 

PUNCH  13,  (J,  (P(I,J) ,I=Ml,M2) , J=l,35) 

Ml=M2+l 

15  M2=M2+7 
GO  TO  1 
END 
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Summary 


Deep  underground  blast-resistant  cylindrical  structures  can 
be  most  easily  designed  by  considering  static  loads  which  are 
adjusted  by  a  factor  to  compensate  for  dynamic  effects.  The 
analysis  presented  permits  one  to  compute  stresses  and  radial 
displacements  in  thick  cylindrical  shells  resulting  from  any 
exterior  static  loading. 
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CHAPTER  2 


DYNAMIC  LOADING  OF 
THICK  CIRCULAR  ARCHES  AND  SHELLS 


Kinetic  Energy  Expression 


The  expression  for  the  kinetic  energy  of  an  arch  Megment 
deforming  in  its  plane  is  readily  obtained  by  integration.  De¬ 
noting  the  kinetic  energy  by  E*  , 


*  -ts!  « «. 


0  A 


in  which  a  is  the  subtended  angle,  p  is  the  mass  per  unit  volume, 
and  w  and  v  are  now  functions  of  time  as  well  as  of  z. 

Substituting  for  v(z,t)  from  equation  (1-1),  letting  w(z,t)  = 
w(t)  and  integrating  with  respect  to  z  over  the  area, 


In  this  expression  I  is  the  moment  of  inertia  of  the  cross 
section  with  respect  to  the  x  axis  (Fig.  1)  ,  and  J  is  the  third 
moment  of  the  area  with  respect  to  the  same  axis;  i.e., 


J  will  be  zero  for  cross  sewtions  having  two  axes  of  symmetry. 


Free  Vibrations 


The  equations  of  motion  of  a  th:ick  circular  ring  segment  de 
forming  in  its  plane  may  be  obtained  readily  from  Hamilton's 
principle, 

6^(E*-V)  dt  *  0. 
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With  Ek  from  equation  (1-20)  and  V  from  equation  (1-3)  this  becomes 

•S‘ *$■*<&■]♦« *  &■] 


'0  0 

_  _ j  ?a 

2 a8  Kbt  ~  B0Bt 

With  the  notation 

I'  =  l/aa A  and  J'  =  J/a®  A 


,  (*2L  .  3a. a .  M  fav  a_  AEZ 

'a*.  2a  lB0  '  2a 


(w  +  |rr)a}d8dt  =  0 
de  J  (1-21) 


the  Euler  equations  of  this  integral  are 

(1+3I'+J',|^--{2I'+J')g2?  -^gf  +f?>=°  *1-22* 

and 

alw-j2®i(I^+JMalw _ +  2^  +  i!”  +  W(1  +1)  +  i  in 

ae4  ez  (I  J  ,deadta  aea  ez  at3^  z'  z  ao 


+  uMa*v  =  0 

ez(21  J  'aoat3  °* 


(1-23) 


It  should  be  noted  that  a  solution  of  equations  (1-22)  and 
(1-23)  is 


v  =  0 


w  =  C*  Sin  J (1+Z)^  +  CaCos  (1+Z)— 


(1-24) 


in  which  Cx  and  Ca  are  arbitrary  constants.  This  solution  repre¬ 
sents  a  motion  such  that  the  ring  remains  circular  and  deforms  with 
time  with  a  circular  frequency  of 


Wo 


-  i  ^(1+Z)  E/p. 


Thus  when  —  V  (1+Z)  E/p  increases  by  an  amount  of  2tt,  the  ring 

d 

undergoes  one  cycle  of  this  motion.  The  period  of  this  motion  is 
2rr  2na 


To  = 


Wf 


V(l+z)  E/p  * 


(1-25) 


The  duration  of  a  transient  load  may  be  conveniently  expressed 
in  units  of  T0 . 
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In  what  follows  we  will  be  interested  in  the  complete  ring  and 
in  motion  which  is  symmetrical  at  all  times  to  the  axis  0  =  0 ;  thus 
we  assume  the  following  solution  to  equations  (1-22)  and  (1-23) . 


w  =  Cos  n0  Sin  ^  V E/p 

®  (1-26) 

v  =  r  Sin  n9  Sin  •Je/o 

In  equations  (1-26)  n  must  be  an  integer  to  satisfy  the  periodicity 
requirements  on  v  and  w,  ^  ^E/p  is  the  circular  frequency  u>,  and  r 

cl 

is  the  ratio  of  the  amplitude  of  the  displacement  v  to  that  of  w. 
Substituting  equations  (1-26)  into  (1-22)  and  (1-23) 

[(2l'+  j')npa-n]  +  r[pa (1+31 '+J ' ) -na  ]  =  0  (1-27) 

rz(na-l)a-pa+l  -  (l'+J')pana"|  +  r[’n-(2l/+J/)pan"]  =  0 

J  '  J  (1-28) 

For  equations  (1-26)  to  be  a  solution  tfo*  equations  (1-22)  and  (1-23)  , 
it  is  necessary  that  the  determinant  of  the  homogeneous  equations 
(1-27)  and  (1-28)  vanish.  Setting  this  determinant  equal  to  zero 
yields 

p*  +  (1+3I'+J')  •  i) 

-pa{(l+3l'+j')[(nsl-l)S+  |]  +  n‘-  |(2I'+J'-  |)na} 

a 

+  na  (na-l)  =  0  (1-29) 

Equation  (1-29)  may  be  used  to  compute  the  natural  fre¬ 
quencies  of  vibrations  of  complete  thick  rings  for  each  mode, 
n=0, 1, 2, . . . For  n=0  it  yields 

p2  =  1  +  Z 

the  motion  described  by  equations  (1-24)  .  For  n=l  it  yields  but 
one  value  of  p8  , 

a _ 2 

P  =  1+41 '+2J '-I /a  * 
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For  n*2,  equation  (1-29)  yields  two  values  of  pa  for  each  integer 
n.  The  higher  value  represents  the  dimensionless  circular  fre¬ 
quency  of  a  deformation  which  is  primarily  extensional  and  the 
lesser  value  of  pa  is  that  of  a  deformation  which  is  primarily 
flexural. 


Equation  (1-29)  is  rigorously  consistent  with  the  Winkler 
assumptions;  it  thus  includes  the  effect  of  "rotary  inertia"  but 
not  the  effect  of  shearing  deform* tion.  Barer?  derived  a  similar 
equation  for  the  frequencies  of  thin  rings  which  accounts  for 
the  coupling  between  the  extensional  and  flexural  motions.  For 
thin  rings  Love4  cites  the  formulas 


and 


uu“ 


ur 


AE  /i  2  \ 
=  — a  (1+n2 ) 
ma 


El  na  (na-l) 
ma4  na+l 


for  extensional  and  flexural  vibrations  respectively.  Here  uu,  A, 
E,  I  and  n  denote  the  same  quantities  as  in  this  paper  and  m  is 
the  mass  per  unit  length  of  centerline.  Equation  (1-29)  includes 
that  of  Baron  and  those  cited  by  Love  as  special  cases. 


A  review  and  survey  of  the  implications  of  the  analysis  pre- 
.  .  ted  is  in  order.  Eqv'.tions  (1-22)  and  (1-23)  are  the  equations 
of  motion.  Solutions  of  the  form  of  equations  (1-26)  exist  pro¬ 
vided  equations  (1-27)  and  (1-28)  are  satisfied.  Nontrivial 
solutions  to  equations  (1-27)  and  (1-28)  exist  for  each  integer 
value  of  n  provided  equation  (1-29)  is  satisfied.  For  both  n  =  0 
and  n  *  1  there  is  a  value  of  pa  which  satisfies  equation  (1-29) . 
For  n  >  1  there  are  two  values  of  pa  which  are  roots  of  equation 
(1-29) .  In  what  follows  these  two  values  will  be  designated  as 
pj  and  pj  where,  arbitrarily,  pj  is  the  lesser  and  pj  is  the 
greater.  For  each  value  of  pa  either  of  equations  (1-27)  or 
(1-28)  will  yield  a  value  of  r.  For  pa  *  pj  the  value  of  r  will 
be  designated  as  rn ;  for  pa  =  pj  the  value  of  r  will  be  called  rB  . 
For  computational  purposes  the  simpler  equation  for  rB  will  be 
used. 


r  =  n-(2i'+j')n  p; 

"  pj (l+3l'+J')-na ‘ 


(1-30) 


In  the  following  analysis,  which  is  meant  to  pertain  to  a 
segment  of  a  long  thick  cylindrical  shell,  J'  will  be  set  equal 
to  zero  and  ro  will  be  given  by 


n  (1-21 *pa ) 

"  pa  (1+31 )-na 
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Since  the  differential  equations  of  motion  are  linear,  solu¬ 
tions  may  be  summed  and  the  resultant  normalized  expressions  for 
w  and  v  may  be  written  as 

CD 

w  =  Sinu!0t  +  Cos0  Sinuul  t  +  ^CosnQ  (Sinuu,  t+SinuuB  t) 

2 


v  =  r,  sine 


Sinuux  t 


CD 

l 


Sin  n0  (r.  Sinuu.  t-r_  Siniw.  t) 


by 


The  uu's  in  these  equations  are  related  to  the  p's,  of  course, 


uj. 


V E/e 

Q 


and  uu0  and  uu1  are  written  as  uu0  and  u)x  so  that  consistently  the 
5.  's  are  the  frequencies  of  the  motions  which  are  primarily  ex- 
tensional. 


GENERALIZED  COORDINATES  AND  THE  LAG  RANG  IAN  FUNCTION 


Knowing  the  modes  of  vibration  and  their  associated  frequen¬ 
cies  it  is  now  possible  to  proceed  with  the  problem  of  the  response 
of  the  ring  to  the  time-dependent  loads.  The  generalized  coordi¬ 
nates  will  be  designated,  following  Baron3  ,  by  q,,  (t)  and  qB  (t)  in 
solutions  of  the  form 


w 


00 

Lk (t)  +  qB  (t)  Jcos  n0 
n=0 


CD 


V  =  Y  [r.q.  (t) 


+  r. 


5,  (t>] 


Sin  n9 


(1-31) 


in  which  it  is  understood  that  qo  (t)  =  qx  (t)  *  0. 

Substitutirg  the  expressions  for  w  and  v,  equations  (1-31), 
into  equation  (1-20) ,  letting  a  «  2n,  and  utilizing  the  orthogon¬ 
ality  of  the  sine  cosine  set 
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E*  *  {C1+r?+3l/rJ+4nrr.  I'  +  l'n3]  (|^  ) 


n=0 


+  2[l+r.E,+3r.?.I'+2(r.+?11)nl'  +  I  'na )  (|*  )  <||. ) 


+  [l+?+3I'rf+4nrl'+l'nal  (1-32) 

□  t  J 

Similarly,  from  equation  (1-3), 

v  =  §*T  7  {[l+nsra+2nr,,+Z(l-na)a]<^ 

+  2^1+na  r„  r,  +  (r,  +r,  )n  +  Z(l-na)  ]q,q„ 

+  £ l+na  ra  +2nrn  +Z  ( l-na  )  *  Jq£  ]■  (1-33) 

In  equation  (1-32)  the  coefficient  of  is  identically  zero. 

ot  o  t 

Likewise  in  equation  (1-33)  ,  the  coefficient  of  qnqB  is  zero. 

(The  proof  that  these  quantities  vanish  is  extremely  tedious  and 
is  not  presented).  Thus,  in  final  form, 

E*  _  {[l+r;+3l'rj+4nr,l'+l'na  j(||*)a 

n=0 

+  [l+rJ+3I/r;+4nrBI/+I/na](|^®)  ]>  (1-34) 

00 

v  =  ^2a  ^[1+narJ+2nr.+Z(l-na)  ]qj 

+  ("  l+narJ+2nrB  +  Zd-n5)’]^  (1-35) 


and  L,  the  Lagrangian  function,  equals  E*  -  V. 


l 

1 
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I  o* 


There  remains  the  task  of  determining  the  generalized  forces 
corresponding  to  each  qBand  .  Once  this  is  accomplished,  La¬ 
grange's  equations  may  be  used  to  deduce  the  differential  equations 
vdiich  describe  each  qo  and  q^  . 


GENERALIZED  FORCES 

The  ring  shown  in  Fig.  6  below  is  considered. 


Figure  6.  Ring  with  time-dependent  loads. 


It  is  loaded  at  its  outer  surface  by  a  radial  load,  X(e,t)  and  a 
tangential  load  Y(0,t). 

To  determine  the  generalized  forces  corresponding  to  q,,  and 
n  ,  X  and  Y  will  first  be  expanded  in  even  and  odd  Fourier  series, 
i .  e. , 
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X(9,t) 


=  +  y  xB  (t)  cos  ne 

n=l 


Y(9,t)  =  ;  Yn  (t)  Sin  n9 

n=l 


(1-36) 


The  virtual  work  of  the  external  forces  X  and  Y  will  be  desig¬ 
nated  by  6We  and 

2n 

6We  =  \  [x6w  +  Y6v)  IRq  d9  (1-36) 

~  z=c 

When  q,,  is  varied  by  6^  ,  from  the  first  of  equations  (1-31)  , 
6w  =  6qBCos  n0 

From  equation  (1-1)  and  equations  (1-31) 


v)  =  -~  v 

z=c  a  a  38 


6v)  =  -i»6v  -  -  6  {“■) 

z=c  a  a  ae 


=  JBov  .£3W 
a 

=  So 
a 

6v)  =  (f°r,  +  -f^lsin  n96q, 

6We  -  nR<,[x,  (t)  +  Y„  (t)  (|°r„  +  -^S-> q. 

F.  =  "Ro[x.  (t)  +  Yn  (t)  (““t,  +  71)] 
in  which  FB  is  the  generalized  force  corresponding  to  q,,  . 

Similarly  the  generalized  force  corresponding  to  q,,  is 

F.  =  "Bo  {*»  <t)  +  Y.  (t)  (f°  ?.  +  ^£)  } 


Thus 

and 


(1-37) 


(1-38) 


The  Lagrangian  function,  L,  is  now  known  from  equations  (1-34) 
and  (1-35) 

L  =  E*  -  V 


and  the  generalized  forces  from  equations  (1-37)  and  (1-38)  .  The 
Lagrange's  equations  of  motion  are  thus 
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aL  _  d  dL 

aq«”  dt  dq„ 

dL  d  dL 
dq«  dt  aq. 


Application 

Two  significant  problems  have  been  worked6 . 

The  first  problem  considered  was  that  of  the  thick  ring,  initially 
at  rest,  loaded  by  a  uniform  external  pressure  which  varied  with 
time  as  shown  below. 


This  loading  produces,  of  course,  only  radial  motion.  Figure  7 

shows  the  radial  displacement  w  as  a  function  of  dimensionless 

time  t  for  various  values  of  the  dimensionless  decay  time  r  , 

d 

t  is  the  ratio  of  actual  time  to  the  period  of  the  n  *  0  mode  of 
vibration  (equation  (1-25)). 

The  second  loading  chosen  was  such  as  to  excite  the  funda¬ 
mental  bending  mode  of  vibration  (n  *  2) ,  and  produce  large 
bending  stresses.  A  radial  pressure  which  varies  as  Cos  29  and 
linearly  with  time  from  P0  at  t  «  0  to  0  at  t  *  was  chosen. 

The  dynamic  response  of  a  ring  with  a  thickness  to  mean  radius 
ratio  of  one  to  three  to  this  loading  is  given  in  Figures  8,  9, 
and  10. 
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Conclusions 


The  analysis  presented  allows  one  to  predict  the  response  of 
thick  rings  to  dynamic  loads.  The  theory  can  be  used  to  determine 
stresses  and  displacements  in  long  cylindrical  shells  imbedded  in 
acoustic  or  elastic  media. 

Since  experiments  in  dynamic  loading  are  difficult  to  perform 
one  can  judge  a  theory  only  by  comparing  its  predictions  with  the 
results  of  other  theories. 

Experiments  show  that  the  formulas  cited  by  Love  (page  37) , 
the  work  of  Baron3 ,  and  equation  (1-29)  all  predict  the  lower 
flexural  natural  frequencies  of  thick  and  thin  rings  very  ac¬ 
curately.6  Since  the  largest  displacements,  strains  and  stresses 
occur  with  the  lower  modes,  the  work  of  Baron,  though  admittedly 
applicable  to  thin  shells,  should  yield  results  applicable  with 
engineering  accuracy  to  thick  shells. 
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SECTION  II 


STRESSES  AND  DISPLACEMENTS 
IN  A 

CIRCULAR  CYLINDER 
FOR 

DISTRIBUTED  STATIC  LOADS  APPLIED 
PERPENDICULARLY  TO  CYLINDER  AXIS 


INTRODUCTION 

The  basic  theory  of  the  method  of  solution  used  here  is  de¬ 
scribed  by  N.  I.  Muskhelishvili1  ,  though  he  does  not  claim 
originality  for  some  of  the  background  work.  The  procedure  is 
restricted  to  the  solution  of  plane  problems  of  elasticity  and 
for  static  cases.  The  governing  differential  equation  for  the 
plane,  static  elasticity  problem  is  known  as  the  biharmonic 
equation,  v*cp  *  0,  where  cp  is  a  stress  function.  It  can  be 
shown  readily1'  a  that  an  equivalent  solution  for  stresses  and 
displacements  can  be  accomplished  by  making  use  of  two  stress 
functions  of  a  complex  variable.  Ultimately  these  two  stress 
functions  are  written  as  Laurent  series'  for  the  particular 
problem  of  the  circular  cylinder.  Then,  with  the  applied  loads 
expressed  as  complex  Fourier  series',  coefficients  of  like  terms 
are  matched  at  the  two  boundaries,  providing  sufficient  equations 
to  solve  for  all  coefficients. 

Figure  11  illustrates  a  cross  section  of  the  circular  cylin¬ 
der  with  an  indication  of  the  loads  applied.  Both  radial  and 
tangential  loads  may  be  applied  to  either  or  both  boundaries;  to 
avoid  confusion,  only  representative  separated  sections  of  the 
distributed  loads  are  shown.  The  total  loading  must  of  course 
be  self-equilibrating  in  order  to  have  a  static  problem,  but  there 
is  no  necessity  that  the  loads  in  each  individual  boundary  have 
zero  resultant.  In  a  later  section,  a  test  for  force  resultant 
and  another  for  moment  resultant  are  presented;  these  tests  are 
included  in  the  computer  program  to  ensure  that  no  attempt  will 
be  made  to  solve  an  improper  problem  by  this  procedure. 

One  precaution  regarding  loading  should  be  made  here;  con¬ 
centrated  loads  cannot  be  handled  by  the  computer  program  as 
developed;  however,  a  concentrated  load  can  be  approximated  by  a 
statically  equivalent  distributed  load  applied  over  an  arc  as 
small  as  one  subdivision  of  the  ring. 

The  input  to  the  procedure  is  in  the  form  of  shear  and  normal 
stresses  on  each  boundary  at  any  desired  number  (this  number  must 
be  divisible  by  four)  of  equally  spaced  (angular  spacing)  points. 
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The  output  produces  radial  (ar) #  tangential  (a^) «  and  shear  (t  ) 

stresses  as  well  as  radial  (u)  and  tangential  (v)  displacements 
at  any  desired  point  of  the  ring. 


Sign  Convention  for  Stresses 

Figure  11.  Cross  section  of  circular  cylinder  with 
portions  of  distributed  loading  indicated. 
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CHAPTER  1 

THEORETICAL  DEVELOPMENT 


FUNDAMENTAL  EQUATIONS 

Since  it  is  well  known  that  the  expression  of  stresses  and 
displacements  for  the  plane  static  elasticity  problem  can  be 
accomplished  in  terms  of  two  stress  functions  of  a  complex  vari¬ 
able,  we  shall  not  repeat  that  development  but  will  start  with 
such  expressions. 


Stresses 

For  polar  coordinates,  the  three  stress  components  are  in¬ 
volved  in  two  equations. 


where 


+  aQ  =  4  Re  $ (z)  =  2^f (z)  +  F(z)J 

(2-1) 

-  ar  +  2i  TrQ  =  2[5#'(z)  +  Y (z) ]e8 1 ® 

(2-2) 

i  =  V-T 

r,0  =  polar  coordinates  (Fig.  11) 

z  =  x+iy  (complex  variable  in  rectangular  coordi¬ 
nates) 


=  re1 9  (complex  variable  in  polar  coordinates) 
z  =  x-iy  (complex  conjugate  of  z  in  rectangular 
coordinates) 


Re 

Hz)  ,Y(z) 

F  (z) 
$'(z) 


re"1 9  (complex  conjugate  of  z  in  polar  coordi¬ 
nates) 

re* 1  part  of 

undetermined  stress  functions  of  complex  vari¬ 
able  z 

complex  conjugate  of  i (z) 
d§  (z) 
dz 


By  subtracting  Eq.  (2-2)  from  Eq.  (2-1)  we  obtain  a  very  useful 
relation  which  does  not  contain  o0 . 

°r  -  1  Tr0  -  Hz)  +  ?(z)  -  e*‘9[5  #'(z)  +  r(z)]  (2-3) 


Relation  Eq.  (2-3)  is  precisely  the  desired  one  for  expressing 
the  boundary  stresses  because  we  see  that  the  radial  stress  is 
the  real  part  and  the  shear  stress  is  the  negative  of  the  imagi 
nary  part  of  the  right  side  of  Eq.  (2-3) . 
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Displacements 

Again  in  polar  coordinates*  the  displacements  can  be  ex¬ 
pressed  in  complex  form. 

2\±  (u+iv)  =  e”1  cp  (z )  -  z  ?(z)  -  ?(z)J  (2-4) 

where 

u  =  radial  displacement  (+  outward) 

v  =  tangential  displacement  (+  in  direction  of  increasing 

0) 

U  =  one  of  the  Lame'  constants  ■  shear  modulus 
a  *  Poisson's  ratio 

cp  (z)  *  (z)  dz 

✓w  -  ^ 

qftz)  =  complex  conjugate  of  cp'(z) 

♦  (z)  *  (z)  dz 

k  a  the  other  Lame*  constant,  k  =  3-4a  for  plane  stress* 
h  =  -7-^-  for  plane  strain. 

1-CT 

We  see  that,  once  we  have  cp  (z)  and  (r  (z),  the  radial  displacement 
is  the  real  part  and  the  tangential  displacement  is  the  imaginary 
part  of  the  right  side  of  Eq.  (2-4)  after  division  by  2 ji. 


Two  Relations  from  Unique  Displacements 

Since  the  cross  section  of  the  hollow  circular  cylinder  is 
a  multiply  connected  region,  the  requirement  that  displacements 
be  single  valued  upon  traversing  a  closed  contour  containing  the 
inner  boundary  results  in  two  essential  relations.  In  general, 
it  has  been  demonstrated1  for  multiply  connected  regions  that  the 
functions  cp  (z)  and  if  (z)  used  in  Eq.  (2-4)  can  be  expressed  as 

A*  In  (z-zk )  +  )  Yk  ln(z-zk)  +cp*(z)  (a) 

k=l 

Hz)  =  )  Yk  In  (z-zk )  +  t*(z)  (b) 

k=l 

where 

m  *  number  of  interior  holes 
zk  *  any  point  inside  the  kth  hole 
Ak  ■  real  constants 
Yk  »Yk  *  complex  constants 

cp* (z)  ,i*(z)  ■  holomorphic  (single  valued)  functions  of  z 
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If  cp  (z)  and  i|i  (z)  from  Eqs.(a)  and  (b)  are  substituted  into  Eq.  (2-4)  , 
we  obtain 

2u£u+ivJ  =  2nie“*  \k+1) ^  z+KYk +Yk  J  (c) 

Lk 

where  |  u+iv  I  =  the  increase  in  (u+iv)  obtained  during  one  tra- 
Ly'  verse  around  a  closed  curve  containing  the 
hole . 

Since  we  must  have  single  valued  displacements,  the  right  side  of 
Eq.  (c)  must  vanish;  therefore,  it  is  obviously  necessary  and  suf¬ 
ficient  that  the  two  relations  of  Eq.  (d)  must  hold. 

Ay  =  0,  KYy  +  Yy  =  0  (d) 

for  k  =  1 ,  2,  . . . ,  m 

In  our  particular  case  of  the  hollow  circular  cylinder,  m  =  1 
(one  hole)  ,  and  we  obtain  the  equations  (2-5)  . 

A  =  0,  Hy  +  y '  =  0  (2-5) 

where  y'  -  complex  conjugate  of  y* 


Since  *  (*)  ■  and  Y  (z)  »  ,  we  can  obtain  Eqs.  (e)  and  (f) 

for  our  case  (m«=l,  zk=0+i0,  z=re* 9)  from  Eqs.  (a)  and  (b)  . 


$  (z) 


-  A(l+ln  z)  +  X  +  SllZL 

z  dz 

=  Alnz+A+^  +  ^ 


(e) 


=  A  In  z  +  t*(z) 
y  dm* (z) 

in  which  A  +  -*-  +  —  is  single  valued, 

z  dz  3 


f  (z) 

Y  ' 

in  which  + 
z 


l!  d»* (z) 
z  dz 

is  single  valued. 


(f) 


Two  important  observations  should  be  made  from  Eqs.(e)  and  (f)  : 
(1)  $  (z)  and  Y  (z)  are  single  valued  (since  A  s  0) ,  and  (2)  y  and 
y'  are  the  coefficients  of  z"1  in  the  series  expressions  for 
$  (z)  and  Y  (z)  . 

Equation  (2-5)  expresses  the  two  relations  referred  to  in 
the  heading  of  this  paragraph.  Actually,  since  the  second  rela¬ 
tion  is  in  terms  of  complex  constants,  Eq. (2-5)  gives  us  three 
scalar  relations. 
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SERIES  SOLUTION  FOR  8TRESSES  AND  DISPLACEMENTS 


Four  It  Representation  of  Load* 

Referring  to  Fig.  11/  we  will  label  the  inner  boundary  of 

radius  as  Lj  and  the  outer  boundary  of  radius  R2  as  La  .  The 

region  between  Lj  and  La  will  be  called  S.  It  is  assumed  that 

the  loading  for  this  problem  will  be  given  in  the  form  of  values 

of  a  and  t  ,  on  L,  and  L- ,  either  as  functions  or  9  of  as  values 
r  rC  ^  3 

at  discrete  points.  Of  course,  as  mentioned  earlier,  the  ring 
must  be  in  equilibrium  from  these  loads  and  the  wieght  of  the 
ring  itself  is  ignored. 


We  shall  represent  the  stresses  acting  on  Lt  and  La  by  com¬ 
plex  Fourier  series. 


(o 


r”iTr8)L1  =  IA‘e'V6  =  I(*+iC*,e 

*  _ _ m 


1  k  0 


—  00 
CO 


—  00 
00 


(2-6) 


‘V^re'L,  =  =  Z<P»+iv*>e' 


k  9 


where 

Ak  =  hit+iCk  =  complex  Fourier  coefficient  determined  by 

loading  on  inner  boundary. 

Ak  =  Pk+;i-vk  =  complex  Fourier  coefficient  determined  by 

loading  on  outer  boundary. 


Laurent  Series  Representation  of  Complex  Stress  Functions 

Though  a  digital  computer  cannot  work  with  complex  numbers 
other  than  by  separation  into  real  and  imaginary  parts,  it  will 
be  convenient  to  continue  the  development  in  complex  form  and 
separate  at  the  end.  A  later  section  deals  in  detail  with  the 
load  representation  by  Fourier  series,  so  for  now  we  will  con¬ 
sider  the  complex  Fourier  coefficients,  A^  and  A^' ,  as  knowns. 
From  Eqs.  (2-3)  and  (2-6)  we  can  express  the  boundary  conditions 
in  terms  of  stress  functions  $ (z)  and  y (z) . 

$  (z)  +  ?(z)-  e3  *  e£z$'(z)+y  (z)  J  = 


For  reference,  one  form  of  equation  (e) 
f  (*)  ■  A  In  z  +  t*(z) 


00 


eik  0 

-®  on  Lj 


(r=R1 ) 


no 

^Aje** 


9 


on  La 

(r=Ra ) 


is  repeated  here, 
(repeated) 


(2-7) 


(e) 


where 


A  ■  real  constant 
$*(z)  *  holomorphic  function  of  z 

Since  $*(z)  is  holomorphic  in  S  up  to  and  including  the  boundary, 
it  can  be  represented  by  a  Laurent  series  (power  series  in  com¬ 
plex  variable) . 

CD 

$*(z)  =  ^z*  (g) 

—CD 

where 

ak  «  complex  constant 


Let  us  now  consider  Eq.  (2-2)  .  Since  the  stresses  are  given  by 
Fourier  series,  they  are  holomorphic  functions;  thus,  the  right 
side  of  Eq.  (2-2)  must  be  also.  From  Eq.  (e)  we  obtain  the 
following: 


*'<z>  =  i  +  ^2i  =  f  +  Ika*z‘"’ 

— OD 


z$'(z)  =  Ar~  + 
z 


oo 

Jka*  zzk 


A  -Q—  +  ^ka*  re"1  9rk  "l  e1  (k  “l )  9 


—  00 
00 
in 


=  A  e~2 1  ®  +  2Jcak  rk  e1  (k  )  0 


—00 


Thus,  we  see  that  z$'(z)  is  holomorphic  and  for  that  reason  so  is 
Y  (z)  .  Therefore,  Y(z)  also  can  be  represented  by  a  Laurent  series. 


w 

$  (z)  =  A  In  z  +  ^a*  z* 

—00 

00 

Y  (z)  =  ^ak'zk 


(2-8) 


Recursion  Relations 

Recalling  Eqs.(2-5)  and  noting  that  y  and  y'  are  the  co¬ 
efficients  of  z"1  in  the  series  representation  of  i  (z)  and  Y (z) 
respectively,  we  see  that  a_x  and  aix  play  the  same  roles  in 
Eqs.  (2-8)  as  do  y  and  y'  in  Eqs.(2-5).  We  can  write  Eqs.  (2-9) 
immediately  from  Eqs. (2-5). 

A  =  0,  na.j  +  a^  *  0  (2-9) 

where 

a^  =  complex  conjugate  of  a^t 
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We  now  substitute  the  expressions  for  *  (z)  and  y (z)  from 
Eqs. (2-8)  into  (2-7) ,  using  A  *  0,  to  obtain  Eqs.  (2-10) ,  two 

very  important  relations. 

00  00  00 

£<l-k)  a,  r*  e1  *  9  +  Ts„  r‘ e-'"  9  -  Ta,'.,r1'-Se'>9 

—eo  woo  — ® 


for  r 


for  r 


Rg 


(2-10) 


Since  Eqs.  (2-10)  must  hold  on  1^  and  La  for  any  value  of  O*0£2tt, 
we  may  equate  the  coefficients  of  like  power  11  of  el0  for  r  * 
and  for  r  »  Rg  to  obtain  recursion  relations  for  the  unknown  com¬ 
plex  coefficients  a*  and  a^  in  terms  of  the  known  complex  coeffi¬ 
cients  A^  and  A£  (remember  that  A*'  and  A£  are  the  complex  Fourier 
coefficients  representing  the  loading) . 


Comparing  terms  independent  of  0  (k=0) ,  we  obtain 
ao  +  io  -  aiaRfa  =  Aj 
from  the  first  of  Eqs.  (2-10)  ,  and 

*0  +  "  a-a  Rg*  =  *o 

from  the  second  of  Eqs.  (2-10).  Since  ag  and  Sg  are  complex  con¬ 
jugates,  these  two  relations  may  be  expressed  as  in  Eqs.  (2-11). 

2  Re  Sq  -  a^gR^*  *  Ag'  (2-11) 

2  Re  ag  -  aiaRJa  =  A£ 

Solving  Eqs.  (2-11)  simultaneously,  we  obtain  Eqs.  (2-12) . 


Re  Bq 


^*0*  -  RjAp' 
2<R?-R?) 

R?«S  (Ab'-Aq') 
Rj-R? 


(2-12) 


Let  us  now  assign  symbols  for  the  real  and  imaginary  parts  of  the 
various  coefficients  involved  here: 


*»' 

K 


ak  +iSk 
Yk+i6k 

hk+iCk 

pk+ivk 


where  k 


0,  - 


1, 
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Thus,  from  the  first  of  Eqs.  (2-12)  we  write  Eq.  (h)  which  will  be 
used  later  in  making  a  check  for  equilibrium. 

Im(RfA3'  -  RfA^)  =  Ra v-  -  R?Co  *  0  (h) 

At  this  point  we  may  wonder  about  the  contribution  of  the  imagi¬ 
nary  part  of  an ;  that  is,  B0 .  The  addition  of  an  imaginary 
constant  to  *  (z)  in  Eqs.  (2-1)  and  (2-2)  leaves  the  stresses  un¬ 
changed;  therefore,  p0  can  be  assigned  any  desired  value,  say 
zero . 

Equating  coefficients  of  e! k ®  for  k  =  -  1,  -  2,  +  ...  -  ®, 
on  r  ~  R.  and  on  r  =  R2 ,  we  obtain  Eqs.  (2-13) . 


( 1-k)  ak  R( 

+  a-.RT' 

-  a,'_2Rr*  =  A.; 

(l-k)a..R- 

-  =  A.; 

We  eliminate  ak'_~  between  Eqs.  (2-13)  to  obtain  Eq.  (2-14)  for 
k  =  -  1,  -  2,  +  ". ..  -  ® . 


(1-k)  (R|-R?  )  ak  +  (RJSk+2-Rr2k+2)a_y  =  A^'R^ +i?-Ak'Rf k +‘ 


(2-14) 

Now  Eq.  (2-14)  is  of  the  form  A  +  iB  =  C  +  iD  since  ,  a..,  ,  A^'  , 
and  Ak  are  all  complex;  therefore,  A  =  C  and  B  =  D.  Thus  we  can 
write  A  -  iB  =  C  -  iD;  that  is,  we  obtain  a  valid  equation  by 
going  to  the  complex  conjugate  form  of  Eq.  (2-14).  However,  in 
doing  so,  we  obtain  a  relation  in  av  and  a_k  instead  of  one  in 
ak  and  a_k .  Thus  we  shall  replace  k  by  -k  in  the  conjugate  of 
Eq.  (2-14)  in  order  to  obtain  another  equation  in  the  same  two 
unknowns . 

(R|k+5-R^k+2)ak  +  (1+k)  (R|-R?  )  a_k  =  Alk  Ri, +s -A^k  R*  +p 

(2-15) 

Equations  (2-14)  and  (2-15)  can  be  solved  simultaneously  for  ak 
and  a_k  provided  the  determinant  of  the  coefficient  matrix  ^  0. 
We  need  consider  Eqs.  (2-14)  and  (2-15)  only  for  k  =  +  1,  +2,  +3, 
...  since  for  each  of  k  =  -1,  -2,  ...  we  obtain  a  pair  of  equa¬ 
tions  conjugate  to  each  pair  fcr  positive  k  and  thus  obtain  no 
new  information  (if  we  know  ak  ,  we  know  ak ) .  The  determinant  of 
the  coefficient  matrix  is  given  as  Eq.  (2-16) . 
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(1-k)  (R>-Rf ) 


D» 


(Rja»+a 


.Rj-a.+a 


) 


(R*'-+».R»‘+»)  (1+lc)  (R»-I^) 


«  (1-k* )  (Rj-Rf  )*  -  (R§>+S-R?l,+S)  (RJak+a-RJ’,ll+s) 

(2-16) 

Prom  Eq.  (2-16)  we  see  that  Dk  vanishes  for  k  *  0,  -  1 ;  therefore, 
Eqs.  (2-14)  and  (2-15)  must  be  solved  specifically  for  these 
values  of  k.  We  already  have  by  the  first  of  Eq.  (2-12) ,  so 
we  will  need  consider  Eqs.  (2-14)  and  (2-15)  specially  only  for 
k  *  1.  Equations  (2-17)  result  from  Eqs.  (2-14)  and  (2-15) 
respectively  for  k  ■  1. 


A^Ra  -  Aj'R*  *  0 

(R|-Rj )  ai  +  2(R»-R?)aLl  «=  A^Rj-A^R* 


(2-17) 


We  now  solve  for  at  and  ait  by  using  the  second  of  Eqs.  (2-9) 


xa.!  +  a^t  *  0 

and  the  first  of  Eqs.  (2-13)  for  k  ■  1. 
a.jR^1  -  aix  Rf1  =  A*' 


(j) 

(k) 


In  order  to  obtain  identical  unknowns,  we  take  the  conjugate  of 
Eq.  (k)  to  obtain  Eq.  (1)  . 

a-i  -  ali  =  A/R*  (1) 


We  solve  Eqs.  (j)  and  (1)  simultaneously  for  a_x  and  aix  after 
which  we  take  the  conjugate  of  a^t . 


a-i 


_  A,% 

K  +  l 


(2-18) 


(2-19) 


Since  a_x  =  from  Eq.  (2-18)  ,  we  can  now  solve  for  ax  from 

the  second  of  Eq.  (2-17)  . 
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(2-20) 


«» 


Rj-Rl 


2  A,'R, 

(R^+RJ)  (k+1) 


It  can  be  shovn  that  Dk  ^  0  for  |k|  ^  2;  and  thus,  for  k  *  2, 
3,  ...  ®,  we  solve  Eqs.  (2-14)  and  (2-15)  simultaneously  for  ak 
and  a.k  .  Also  notice  that  Dk  *  D_k  . 


(1+k)  (Rj-R?)  (Ak'Rjk+a-Ak,RJk+a)-(R£ak+a-RTak+a)  (A^  Rj  +a -Alk  R$  +a  ) 


(2-21) 


(1-k)  (Rj-R?)  (A:kRk^a-AlkRk-ha)-(Rak^a-Rak-Hi)(A:ic  R^k  +a  Rf*  +a ) 


(2-21*) 


Thus,  Eqs.  (2-18) ,  (2-20) ,  (2-21) ,  and  the  first  of  Eqs.  (2-12) 
completely  determine  all  ak ,  the  coefficients  for  $*(z). 

We  now  shall  determine  a^,  the  coefficients  of  the  series  for 
¥ (z) .  Knowing  ak  and  a.k ,  either  of  Eqs.  (2-13)  may  be  employed 
to  solve  for  a^.g .  Using  the  first  of  Eqs.  (2-13),  we  obtain 
Eq.  (2- 2 2a)  which  can  be  re- indexed  and  expressed  as  Eq.  (2- 2 2b) . 

a;.,  =  (l-k)Raav  +  a_k  R^3k +=-  A„'R^k+a  (2-22a) 

for  k  «  -  1 ,  -  2 ,  ... 

-  -(l+k)R?a,+a  +  a.  (k+a)Rr  (»»+»)  -  \'+aRT‘  (2-22b) 

for  k  *  0,  -  1,  +2,  -  3,  +  ... 

Neither  Eq.  (2-22a)  nor  Eq.  (2-22b)  gives  the  term  a!a ,  but 
we  already  have  this  from  the  second  of  Eqs.  (2-12) ;  therefore  all 
coefficients  for  ¥ (z)  can  be  computed.  In  order  to  reduce  the 
number  of  compu  tat  ions,  it  was  decided  to  let  R,  *  1  and  R1  be  the 
appropriate  fraction  less  them  one;  then  the  solution  for  a  cylin¬ 
der  with  Rg  /  1  (but  geometrically  similar)  can  bo  readily  found 
by  dimensional  analysis. 

The  use  of  Eq.  (2-22a)  or  Eq.  (2-22b)  for  large  k  resulted  in 
lack  of  precision  on  a^  because  of  small  differences  of  large  num¬ 
bers.  It  can  be  seen  that,  since  Rl  <  1,  we  have  numbers  less 
than  unity  to  large  negative  powers.  Because  of  this  precision 
trouble,  the  second  of  Eqs.  (2-13)  instead  of  the  first  was  used 
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to  compute  a*'  since  it  contains  only  R,  ■  1,  and  therefore  no 
difficulties  caused  by  small  differences  of  large  numbers.  Equa¬ 
tion  (2-22c)  also  gives  a*' ,  but  from  the  second  of  Eqs.  (2-13). 

a,'  -  -U+k)a4+,R|  +  a_(t+a)H;(“+a)  -  V+SRT 

=  “  +  —  (k  +s  )  “  +a  (2— 22c) 

for  k  **  0,  i  1,  +  2,  !  3,  +  ... 


Summary 

Coefficients  of  i (z) : 

*o . 

*i . 

a-i . 

ak  (for  k  =  i  2,  -  3,  +  ...)  .... 


Equation  (2-12) 
••  (2-20) 

H  (2-18) 

••  (2-21) 


Coefficients  of  f (z)  : 

. Equation  (2-19)  or  (2- 22c) 

ala .  "  (2-12) 

ak  (for  k  *  0,  ±  1,  +  2,  ±  3,  ...)  .  "  (2-22c) 


Equilibrium  Verification 

Inasmuch  as  the  entire  procedure  here  is  restricted  to  the 
static  case,  there  are  three  relations  which  must  be  satisfied  if 
the  hollow  cylinder  is  to  be  in  equilibrium.  With  the  entire 
ring  as  a  free  body  in  Fig.  12,  we  will  sum  moments  about  the 
origin  and  also  sum  forces  horizontally  and  vertically. 
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Equilibrium  of  Entire  Hollow  Cylinder 


From  Eqs.  (2-6)  ,  we  may  write  the  following  relations  after  expand¬ 
ing  and  A*  into  real  and  imaginary  parts: 

00 

(a  )  -  i  (t  )  =  y(rbt+i  )  (Cos  k0+i  Sin  k0) 

r  Lx  re  Lx  Jaf 


(a  )  -  i  (t  ) 

1  L,  r0  L, 


w 

=  £(p„+i  v* )  (Cos  k0+i  Sin  k0) 


Thus,  we  can  separate  into  real  and  imaginary  parts  and  express 
or  and  Trg  on  each  boundary. 
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'v,  -i 


(i\  Co*  k0  -  Ck  Sin  k0) 


(Tr0)  .£<Ck  Cos  k0  +  Sin  k0 ) 

I»i  -® 


(Of)  *  CO®  ^0  "  vk  Sin  ^0) 

Lg  -® 


(t rg )  ”“7 (vk  Cos  kfl  +  0k  Sin  k0) 


Lg  -« 


We  now  substitute  (t  )  and  (t  )  into  the  moment  equilibrium 
equation.  r0  Lj  r0  Lg 


CD  * 

^Mo  *  -Rg  7  [k*sin  k0“  fc*Cos  k0]  +  ^  Z  [k  Sin  k0~  k* COS  ke] 

mmCO  0  _0C 


2  n 
0 


This  expression  for  M  identically  vanishes  term  by  term  except  for 
k  =  0  which  must  be  handled  separately. 

£m„  -  -H^v.  d6+  R?^nC0  <30  «  2tt(Rj  Co  -  Rj  v„)  -  0  (2-23a) 


•  •  Ri  Co  “  Rg  vo 


(repeated) 


(h) 


We  see  that  the  forced  condition  of  Eq.  (h)  is  a  consequence 
of  the  fact  that  the  entire  ring  must  be  in  moment  equilibrium. 

Similarly,  the  entire  ring  must  be  in  force  equilibrium 

V  F  =  0  =  \2nr  (o  )  Rg  -  (a  )  Rilci  Cos  0  +  j  Sin  0)d0 
J0LrLg  rLi 


+  ^2n[<Tr0J  *8  "  (Tr0)  sin  0+1  Cos  6) 


d0 


Again  we  substitute  the  series  expressions  for  stresses. 


00 


Jp  *  y  n{  7  [  (PkRa-^kRi  )C°s  *9  -  (Vk  Rg  -Ck  Ri )  Sin  ke] 

0  —  OD 

(l  Cos  9+5  Sin  0jd0 
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(vkRa-CkRi  )Cos  k9 


+ 


(pk  Ra-i\  )  Sin  ke  j 


(-1  Sin  0  +  j  Cos  0) 


Now,  since  Sin  n0  and  Cos  nd  comprise  a  set  of  orthogonal  functions 
over  the  range  from  0  to  2tt,  we  have  only  the  terms  for  k  *  1  and 
of  like  kind  survive. 

If  =  2TT(p1Ra-Ti1Rl)i  -  2tt(v1R8-CiRi)  j  s  o  (2-23b) 

pxRa  -  T]x  Rx  =  0,  Vi  Ra  -  CiRi  -0  (repeated)  ( (2-17)  first) 

Thus,  we  see  that  the  two  forced  conditions  (one  from  real  equality, 
one  from  imaginary  equality)  of  Eq.  (2-17)  are  consequences  of  the 
necessary  force  equilibrium  of  the  entire  ring. 

In  the  computer  program  for  the  solution  of  this  problem, 
these  three  verifications  are  made  early  in  the  computation  as 
soon  as  Co  <  v0 ,  px ,  ,  Vj ,  and  Cx  are  available.  Of  course, 

with  loading  data  supplied  numerically,  we  would  not  expect  an 
exact  clieck;  but  to  some  tolerance,  these  three  relations  must  be 
satisfied  before  we  have  a  valid  problem  for  the  method. 


Stresses  in  Terms  of  a,  B,  v.  and  6 

Since  a  digital  computer  cannot  function  d:  rectly  in  terms 
of  complex  numbers,  we  will  now  express  the  three  stress  components 
in  terms  of  real  numbers. 

From  Eqs. (2-1)  and  (2-8) ,  we  obtain  Eq. (2-24)  after  using  A=0 
in  Eq. (2-8) . 


a  +aft 
r  0 


CD  CD 

=  (z)  +?  (z)  ]  ■  2[  £  ak  zk  +  ] 

^  —00  —00 

00 

=  2  ^rk  (akelk9  +  ake“lk®) 

—00 

00 

=  2  Jrk£(ak+i0k)  (Cos  k0+i  Sin  k0)  +  (ak -i0k )  (Cos  k0-i  Sin  k9)J 

(2-24) 


—  00 
ao 


=  4  ^Trk  (akCos  ke  -  0kSin  k0) 
—00 
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We  obtain  equation  (2-25)  by  uaing  Eqa. (2-2)  and  (2-8) . 


o-o  +  2it  . 
0  r  re 


6 


2[z  *  '  (z)  +  T  (*)]••* 

OP  OP 

2e* 1  e[i  ^ica,,  zk  “x  +  ^a^ zk  J 

•08  « <n 

OP  OP 

2e810^re"1®  ^kak rk “l e1  (*-l)9+  ^ak#rketk®J 

•  08  m  QP 

®  00 

2C  J)cakelk®  ♦!<*».«  *-+•>«] 

•  00  —00 

00 

2  ^rk£kakCos  ke-k0kSin  lc0-#-yk  Cos  (lce-*-2e )  —  6k  Sin  (Ic9-i>20 )  J 


—00 

00 


+  2i£rk£k0kCos  ke+kakSin  ke+6kCos  (k0+20)  +yk  Sin  (ke+20)  J 

(2-25) 


—  00 


By  adding  Eq. (2-24)  and  the  real  part  of  Eq. (2-25)  ,  we  obtain  a 
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as  given  by  Eq. (2-26) ;  by  subtracting  the  real  part  of  Eq. (2-25) 
from  Eq.  (2-24)  ,  we  obtain  or  as  given  by  Eq.  (2-27)  ;  and  the  imagin¬ 


ary  part  of  Eq. (2-25)  gives  TrQ  directly  in  Eq. (2-28) . 

00 

2a0  *  £  ^4rk  (akCos  k0-0kSin  k0 )  +2rk £kak Cos  k0-k0k Sin 


ke 


or 


or 


+  YkCos  (k9+29  )-6k  Sin  (k0+20)  J  J- 
00 

o0  *  £rk£(2+k)  (akCos  k0-0kSin  k0)  +  ykCos(k0+20) -6k  Sin(k0+20)  j 
—00 

.  (2-26) 

2or  *  £rk£(4-2k)  (akCos  k9-0kSin  k0) -2vkCos  (k0+20) 

—00 

+  26k  Sin  (k0+20)  J 
00 

or  =  £rk£(2-k)  (akCos  k0-BkSin  k6)-YkCos(ke+20H-6k  Sin(k0+2e)  J 

(2-27) 
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ao 

t  =  ^rk£kpkCos  k0+kakSin  k8+6k  Cos  (k0+20)+Yk  Sin  (k0+29)  J 

(2-28) 


Displacements  in  Terms  of  a,  g,  v,  and  6 

We  find  displacements  u  and  v  from  the  real  and  imaginary 
parts  respectively  of  Eq.  (2-4) . 

2^i  (u+iv)  =  e”1  9£*cp  (z) -zq> '  (z) (z)  (repeated)  (2-4) 

where 

0° 

cp  (z)  =  ^t(z)  dz  =  ^$*(z)+A  In  zjdz  ■  z*  dz 

—  CD 


from  Eq. (2-8)  with  A=0 


-2 


*  Y  zk+l+  a-i  In  z  + 


L  k+1 

—  00 


z‘ 


+1 


00 


cp'(z)  =  I  (z)  =  A  In  z  +  £akzk  =  ^zk  ,  since  A  =  0 

—  00  —00 


00 


<p'(z)  =  ^akzk  =  Ja*  rk  e~!  k  9 


—00  —00 


ZCp 


00 

'  (z)  =  ^rk+le-*  (*-*  )9 


*  (z) 


00 

=  (z)  dz  =  ^  ^ak'zk  dz  from  Eq.  (2-8) 

•00 

=  Ik?t2“+l  +  ln  z*l^z'+l 


ijf(z) 


=  fe**+l  +  5->  ln  5  +  Imz'*' 


The  above  integrations  were  carried  out  without  including  any 
additive  constants  because  such  constants  represent  rigid-body 
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motion.  That  this  is  so  can  be  seen  by  observing  that  the  total 
of  such  constants  in  the  square  brackets  of  Eq. (2-4)  would  be  of 
the  following  form: 

(u+iv)  =  -  -  -  j~a  ;-f  ib;1  =  e"‘  ®  (a+ib) 

due  to  2M  L  J 

additive 
constants 
in  Eq.  (2-4) 

=  (a+ib) (Cos  0+i  Sin  9) 

=  a  Cos  9  +  b  Sin  9  +  i (b  Cos  9  -  a  Sin  9) 
u  =  a  Cos  9  +  b  Sin  9 
v  =  -a  Sin  9  +  b  Cos  9 

As  seen  from  Fig.  13,  these  values  of  u  and  v  are  exactly  the  radial 
and  transverse  components  of  displacement  each  particle  would  have 
if  the  entire  body  were  rigidly  translated  by  an  amount  a  in  the 
x-direction  and  an  amount  b  in  the  y-direction. 


Figure  13.  Radial  and  transverse  displacements  of  a  particle 
due  to  given  Cartesian  displacements. 
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Similarly,  we  find  by  Eqs.  (2-1)  and  (2-2)  that  an  addition  of 
an  imaginary  constant  (i0o)  to  $  (z)  leaves  the  stresses  unchanged. 
This  addition  of  i0o  to  *  (z)  ,  however,  changes  cp  (z)  by  the  amount 
of  i0oz  and  changes  the  displacements  as  follows: 

(u+iv)  =  i  —  k 0O  z  =  ie"1  ®  re1® 

due  to  adding  M  M 

ie0  to  §  (z) 

=  A0  ir 


u  =  0,  v  =  A0r 

where  A0  =  or  0O  = 

2n  k 

However,  the  displacements  u  =  0  and  v  *  rA0  are  exactly  those  each 
particle  of  a  body  would  receive  if  the  body  underwent  rigid  body 
rotation  by  the  angle  A 9  about  an  axis  through  the  origin  as  indicated 
in  Fig.  14.  Actually,  of  course,  the  displacement  v  =  rA0  in  the 
instantaneous  v-direction  is  valid  only  for  small  A0.  However,  the 
path  of  a  particle  is  a  circular  arc  of  length  rA0  even  for  large  A9 . 
Thus,  the  addition  of  the  term  i0o  to  the  series  for  i (z)  has  the 
effect  of  a  rigid  body  rotation  of  the  ring  about  an  axis  through 

the  origin  by  an  amount  A0  =  • 


Figure  14.  Radial  and  transverse  displacements  due  to 

rigid  body  rotation. 

In  programming  this  procedure,  these  rigid  body  displacements 
were  left  unspecified  and  thus  the  computed  displacements  are 
"floating."  This  actually  enables  one  to  obtain  the  displacements 
in  an  advantageous  form  especially  for  problems  of  symmetrical 
loading . 
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Ignoring  these  constants  of  integration  in  Bq. (2-4)  and  the 
term  i0o  in  *  (z) ,  we  obtain  a  series  expression  for  displacements 
as  given  by  Eq. (2-4'). 


2m  (u+iv) 


zk+l  +  In  z  + 


zk  +l 


00 

—  **  e"1  )  9 

—  0D 


In  z 


00  — 


-Z£r*+‘] 


2m (u+iv)  *  e-.9['Jx  rk  +l  e1  (“  +l )  0  +  Ka.l 

-00 


(In  r+ie) 


+  rk  +l  e1  (k  +' )  9  -  rk  +l  e”1  (k  "' )  9 


—  00 


-Jfct rr 


-k+ie-i  (k+i)0  _  (in  r-i0) 


or 


00  — 


-  Z&rr‘+le" (,+l)9] 


2u  (u+iv)  *  )k  rk+lelk®  +  e“*9ln  r(na_l-al1) 


Jk  r‘  +l  e*  *  9  -  Ya,  r‘  +l  e*' « 9 
rt  -® 


-  £&r r 


.k  +1  e- 1  (k  +•  )  0 


—00 


ife" 


+1  p-I  (k  +1  )  0 


+  i0e-‘®(Ka.l+aIl)  (2-4') 

where  (Ha.l+aix)  *  0  by  Eq.  (2-9)  . 

Finally,  we  obtain  the  series  solution  for  u  from  the  real  part  of 
Eq.  (2-4')  and  the  solution  for  v  from  the  imaginary  part. 
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2[iu  -  j[~[ak  Cos  ke-BkSin  kej+ln  (na^ -y_i )  Cos  0  +  (kB_x  -ft.*,  )  Sin  ej 
.00 


+  p£f[a‘co9  k9-e*  Sin  ksj-  pk+l£akCos  k0-0kSin  kQ J 

0  -OB 


k  Cos  (k@+20) -6k  Sin  (k0+20)  J 


j^Yk  Cos  (k0+20 )  -6k  Sin  (k0+20 )  J 


(2-4a) 


-2 


2[iv  =  j— jjik  Sin  k0+0kCos  k@J+ln  rjup^ +6„x )  Cos  0  +  (y-i  -na-l )  Sin  ©J 

.00  ^ 


00  ^  00 
+  ■  ^ak  Sin  k0+f3kCos  koJ+  pk+l£akSin  k0+0kCos  k@J 

0  -« 

+  ^Yk  Sin  (k0+20)  +6k  Cos  (k0+20)  "| 


OD  ^  |  ^ 

+  ^k+1  £Yk  Sin  0*9+20 )  +6*  Cos  (k0+20)  ] 


(2-4b) 


The  form  for  Eqs.  (2-4a)  and  (2-4b)  is  obtained  using  ak  *  ak +ipk 
and  ak  =  Yk  +i&k • 

Since  the  Fortran  system  cannot  handle  negative  subscripts  for 
subscripted  variables,  we  will  later  convert  the  stress  equations 
(2-26,  2-27,  2-28)  and  the  displacement  equations  (2-4a,  2-4b)  into 
forms  not  using  negative  indices. 


Summary 

Stresses 

°e  ’ 

®r  ’ 
Tr8  ' 


Equation  (2-26) 
"  (2-27) 

"  (2-28) 
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Displacements 


u  .  Equation  (2-4a) 

v .  "  (2-4b) 


The  above  equations  for  stress  and  displacement  are  in  terms 
of  series  coefficients  ak  ,  0k  ,  yk ,  and  6k  which  are  the  real  and 
imaginary  parts  of  coefficients  ak  and  ak  .  The  coefficients  ak 
and  ak  are  in  turn  expressed  in  terms  of  A^  and  A£ ,  the  Fourier 
loading  coefficients;  ard  these  loading  coefficients  are  determined 
from  the  distribution  of  a  and  t  .  on  L,  and  La.  Thus,  our  next 

task  is  to  develop  a  useful  procedure  of  determining  Ak'  and  from 
the  given  loading  on  and  La .  The  loading  on  the  boundaries  may 
be  in  graphical  or  numerical  form  and  the  usual  formal  determina¬ 
tion  of  Fourier  coefficients  will  not  suffice. 
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REPRESENTATION  OF  a  and  t  ON  BOUNDARIES 

r  r0 

BY  COMPLEX  FOURIER  SERIES 


Figure  15.  Representative  distribution  of  one  of  the 
stresses  (o  )  on  one  boundary  showing  the 
"double  value"  nomenclature. 


Our  next  task  is  to  spell  out  a  detailed  procedure  for  ob¬ 
taining  the  complex  Fourier  coefficients  A*'  and  A£  in  Eqs.  (2-6) 
which  are  repeated  here  for  convenience. 


GO 


ar  “iTr0  =  I  **  e,k0  on  ^ 


(repeated)  (2-6) 
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where 


A„#  *  i\  +  i£k  ■  complex  Fourier  coefficients  on  inner 

boundary. 

A^'  =  pk  +  ivk  =  complex  Fourier  coefficients  on  outer 

boundary. 


In  the  usual  manner,  we  obtain  the  Fourier  coefficients  by  multi¬ 
plying  both  sides  of  each  of  Eq. (2-6)  by  e”tB0  and  integrating  with 
respect  to  0  from  0  to  2tt. 

2rr  ®  2rr 

\  (o  -iT  )  e“lB9d0  *  V  Ak'  \  e1  ('f“B)0d0  *  2nAk' 


k=-®  "0 


2tt 


lince  ^  e1  (k"")6d0  =  0  if  n  /  k 

=  2n  if  n  =  k 


**'  =  2^„(ViTr9).  e",,0d0  *  *5  "or-iTr0).  (COB  W'iSi"  M) 


as 


i,2"r  -i 

=  -r— \  (a  Cos  k0-T  Sin  k0)  -i(T  Cos  k0+o  Sin  k0)  d0 

l  r  re  T  r0  r  9  J 


2rr 


1  f 

=  —  \  (a  Cos  k0  -  t  ftSin  k0)T  d0 
2n  J  r  r0  I* 


2tt 


(2-29) 


1  r 

Ck  =  -  *5—  \  (t  Cos  k0  +  a  Sin  k0)  de 
J0  0  Lx 


In  similar  fashion,  we  can  write  expressions  for  Ak  ,  pk  ,  and  vk . 


2tt 


K  = 


Pk  = 


—  \  r (a  Cos  k0-r  Sin  k9)  -i  (t  «Cos  ke+a  Sin  ke)  Id© 
2n  L  r  r0  _  r0  r  _  J 

L s  La 


0 

2tt 


-r—  \  (a  Cos  k0-T  „Sin  k0)  d0 
2"  Jp  r  re  ^ 


\  %V- 


2tt 


-  t-  \  (t  Cos  W+  a  Sin  ke)  d0 
2tt  re  r 


(2-30) 
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In  case  the  loading  is  given  in  functional  form,  the  coefficients 
can  be  found  from  Eqs.  (2-29)  and  (2-30)  and  then  read  directly  into 
the  computer  program.  A  subprogram  to  compute  these  coefficients 
could  be:  used.  In  general,  however,  it  is  anticipated  that  the 
loading  stresses  on  the  boundaries  will  be  in  graphical  or  numerical 
form?  the  remainder  of  the  development  of  this  division  will  be 
aimed  at  obtaining  the  coefficients  p,  v,  ri,  Q  for  this  situation. 

A  representative  distribution  of  one  stress  (of)  on  one  boundary 

is  shown  in  5ig.  15.  We  will  represent  the  given  boundary  load 
stresses  as  segments  of  linear  functions  of  6  between  equally  spaced 
divisxons  of  the  boundaries.  These  linear  functions  of  0  would  be 
straight  lines  on  a  rectangular  stress  -  vs.  -  0  plot,  but  do  not 
plot  as  straight  lines  on  one  such  as  Fig.  15.  In  order  to  allow 
for  discontinuities,  each  loao  stress  on  each  boundary  is  assigned 
two  values,  a  "near”  value  and  a  "far"  value.  These  are,  of  course, 
the  values  of  the  considered  stress  just  prior  to  and  just  after 

F 

the  discontinuity  respectively.  In  a  symbol  such  as  (a  )  ?  n  re- 

th  r*  w 

fers  to  the  n  division  point  of  the  boundary,  Lg  refers  to  the 
outside  boundary,  and  F  refers  to  the  value  of  a r  on  the  "far"  side 

of  the  discontinuity.  At  most  division  points,  there  will  be  no 
discontinuity  and  the  "near"  value  of  the  stress  will  be  equal  to 
the  "far"  value. 


Loading  Coefficients  for  Case  where  Discontinuities  May  Exist 


For  each  linear  segment  of  a  loading  stress,  we  may  express 
this  stress  in  the  form  given  by  Eq. (m) . 


_  k„+l>"  -  <vL  K9*’ 

„  .  F  B  L_  *  T i— 

(o  )  + - — - — - 

n  (8.+X-0.) 


(m) 


where 


a  between  0B  and  0B+l  on  outer 
boundary  at  the  angle  0 . 


To  reduce  writing,  we  shall  develop  the  expression  for  pk  only  and 
it  is  understood  that  we  are  working  only  with  loading  stresses  on 
the  outer  boundary.  From  the  first  of  Eqs. (2-30)  ,  we  express  p„ 
by  integrating  over  each  segment,  one  after  the  other,  until  the 
integration  around  the  entire  boundary  is  complete. 
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2tt pk  =  ^  §  *  +  (iar^La  Cos  k9  "  ^Tre)L8  Sin  ke]d0  <n) 

n=0  0B  n  n 

where  M+l  =  number  of  equal  subdivisions  of  the 
circle. 

Into  Eq. (n)  are  now  substituted  the  linear  expressions  for  (o  ) L* 

r  r  n 

and  (t  in  the  form  given  by  Eq.  (m)  .  For  convenience,  we 

r9 

n 

shall  represent  by  a  and  t  by  t;  also  we  understand  that  both 
a  and  t  are  the  loading  stresses  on  Lg . 

k  V9n+l{f0f  (9-en)1cos  k9 

n=0  J9  LL  A9  J 


2no 


N  F 

-  Ft;  +  (T.n.+1.  ■"  T,b  )  (Q-9b  )  Is  in  kejdQ 
L  A9  J  J 

where  9„ +i -9n  =  A9  since  each  boundary  is  divided  into  M+l  equal 

parts. 

M  0.,  N  F  N  F 

=  V  \  Ta  -  en  (a°n  ~  °n )  ~|cos  k9  +  (.?-n±x ■  )  9  Cos  ke 

n^O  1  L  A9  J  A0 


P  N  F 

F  Tn+1  -  T  n 


-  r 

-  _Tn-0n  ( 


A9 


N  F 

Is  in  k9  +  (Ii±i - -J5L)9  Sin  ke]>d0 

J  A  9  J 


M 

n=0 


{  K  -  e„(%^)][^]9”+'+  k9+ ksCos  k<] 


F  .  -T^t  -  T.F,[-C08  (T^,  F 

A9  L  k  J 


♦[xf  -  9.  C 


'n  +i 


A9 


— )[^Cos  kQ-  —g  Sin  ke]  ] 


We  now  evaluate  the  above  relation  for  2npk  in  terms  of  the  indicated 
limits  and  then  cancel  all  possible  terms. 


n  +1 

n 
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2npk  =  7  7  {of+1Sin  k0B+l  -  ofsin  k0B  ~(^^-J-CT|)  (Cos  k0n+l 
n=0  *  A6 

-  Cos  k0B )  +  Tf+l  Cos  k0B+1-  Tf  Cos  k0B 
,  N  F 

-  -(Ll+UILl)  (sin  k0B+1  -  Sin  k0B)|  (p) 

K  A0  J 

for  k  =  i  1,  i  2,  .  .  .  -® 


By  writing  out  several  terms  of  the  above  series  we  see  that  we 
can  reassemble  the  terms  into  a  more  useful  form  as  given  by  Eq. 
( 2— 3 lp )  .  Notice  that  o0  =  oM  +x  arid  t0  =  t.+x. 


Pk 


N 

ar  + 

Sin  k0.  +  ( — s — 


F 

o 

r. 


r«  -i 


-  a 


N 

r 

rn+l 


kA0 


)  Cos  kQ 


J 


2kn 


M+i  „ 

u\- 

n=l  n 


Trft  }  COS  k0n-<- 

r0,I* 


N  F 

T  +  T 


F  N 

“  T  -  T 

r0.  - 


kA0 


- rla+i.)Sin  keB  ] 


for  k  =  -l,  -2 , 


+ 

.00 


(2-31p) 


Since  this  procedure  is  designed  for  use  with  Fortran  programming 
of  digital  computers,  we  must  avoid  the  use  of  subscripted  variables 
with  negative  subscripts.  Therefore,  we  shall  relabel  pk  ,  vk  ,  r|k  , 
and  Ck  in  Eqs. (2-29)  and  (2-30)  for  negative  k  as  £k  ,  ,  and 

Ck .  If  -k  is  substituted  into  Eq. (2-31p)  for  k,  we  obtain  no  change 

in  the  first  ^  but  we  obtain  a  sign  change  for  each  part  of  the 

second  Therefore,  we  do  not  need  to  compute  pk  separately;  if 

we  represent  pk  as 

r'P  vP 


Pk 


Lo  +  U'  (k  *  i,  2'  •••  •) 


then  we  will  find  pk  by 

■P  r'P 


*  -Z-Z-  (k  = 


1,  2, 


(2-31p) 


We  can  compute  vk  and  in  a  very  similar  manner  to  that 
above  for  ok  and  pk  to  obtain  Eqs. (2-31v)  and  (2-31vJ 
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Sin  k0B  + 

L, 


N  F  P 

T*  +  T  —  T 

(  r8,  reB  gfta-i 

kA9 


tn 

-£ii±L)CQ8  keB 
La 


oF  )  Cos  ke, 
r 

1  La 


kA0 


N 

°r 

J JL±i.)  Sin 

L, 


(2-31v) 


(k  -  1,  2,  ...  - 


(2-31v) 


The  coefficients  r|,  -q ,  Q,  and  Q  are  confuted  in  exactly  the 
same  way  except  that  they” are  the  Fourier  coefficients  for  the 
inner  boundary  and  thus  are  computed  from  the  shear  and  normal 
loading  stresses  on  Lx . 


N 

ar 

)  Sin  k0B  +  (-li 
Li 


Hi 


Cos 


N  F  _  F  N 

TP  )  Cos  keB-(-Eis — iii — E5*±L)  Sin 

reBT  “  XA9 ; 

Li  h, 


(k  *  1 ,  2 ,  ...  « ) 


(2-31r|) 


H ** 


(2-31ri) 


N 


N 


re. 


k^e 


k  = 


+ 


F 

o  )  Cos  ke. 
r  “ 

“L, 


kA9 


N 

°r 

— 1±L)  Sin 
Li 


(2-310 
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y-t-t 


(k  «  1,  2,  ...  ®) 


(2-310 


Notice  that,  if  discontinuities  are  present  in  the  loading 

stresses,  the  coefficients  diminish  as  — ■  .  This  is  a  well  known 

k 

result  in  Fourier  analysis. 


Loading  Coefficients  for  Case  of  No  Discontinuities 


At  this  point,  it  should  be  noted  that  the  first  version  of 
the  Fortran  program  was  based  on  the  assumption  that  there  were 
no  discontinuities  in  the  loading  stresses;  for  this  case,  the 
"near"  and  "far"  values  of  any  loading  stress  are  the  same  at 
each  division  point  and  we  can  eliminate  the  superscripts  N  and 
F.  Thus,  for  no  loading  stress  discontinuities  we  obtain  from 
Eqs.  (2-31)  the  following  expressions  for  p,  p,  v,  v,  r),  r|,  Q, 
and  Q  for  k  =  1,  2,  3,  .  . .  ® : 


T 


2nA9ka  L 
M±1 


>  (2t 


2rrA0ka  L>  '  'r9_ 
n=  l 


a  -  o  )  Cos  k0,, 
r»-l  L- 


t  -  t  )  Sin  kQ 

re.-i  re>+l  ^ 


(2-31p)  # 


Y'p  y'p 
^  o  ^  T 


'0  r'o 


(2-310)  ' 


v,  =  - 


1 

L  (2t~  '  T~  '  T 


n=l  r9”  r9-+> 


)  Cos  k8. 


(2'31v)' 


c-'v  r-'V 

L  *  L  „ 
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(2-31v)  ' 


H±1 


TV  *  _  ^  ,  a  )  (2o  -  0  -  0  )  Co*  k0„ 

2ni0k  n4i  r.  r.  -i  r.  +i  V 


_1 _ 

2ni9ka  nZi1'2Tr9,”  Tre,.1 


Tro  >  sin  k«. 
«.+i  ^ 


(2-31t,)  ' 


•'tl  _  y'n 

’  a  “t 


tw 

—  L-*  Q  -J  T 


(2-31r))  ' 


c„ 


1  To 

nil  Tr0.  Trfl*->  Tr0.+> 


2nA0ka 


)  Cos  k0B 


i  77i X  ^  (2o„  **  a*.  “  a  )  Sin  ke 

2nA0k  n-l  r>  r»-i  r»  +i  j. 


(2-310  ' 


-x>i 


'C 

a 


1.-IM 


'C 

t  “a 


(2-310  ' 


Notice  that,  for  no  loading  discontinuities,  the  coefficients 
diminish  as  this  fact  is  also  well  known  in  regard  to  Fourier 
series. 
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Loading  Coefficients  for  k  =  0 


It  is  readily  seen  from  Eq. (p)  that  Eqs. (2-31)  and  (2-31) ' 
are  invalid  for  k  *  0  since  k  appears  in  the  denominators;  this 
situation  was  caused  by  the  integration  of  Eq. (n)  with  t\e 
stresses  in  the  form  given  by  Eq. (m) .  Therefore,  we  must  handle 
the  zero  terms  separately  beginning  with  expressions  in  the  form 
of  Eq. (n)  for  k  *  0. 


or 


2np0  =  ^  \  (ar)^8de 


Similarly, 


n“cre 


M  0..  „  N  F 

M  k  +  -A)(9-en)]  de 


N 


AQ 

F 


=  ^  \on  9  +  (^iii - 2±)  4-  -  9b  9 )  ]  1 

u  L  A9  2  Jfl 

n=0 


n=0 


■l 


{of  (9„  +1  -9, !  +  °M  [b;+>2~  b;  -0.(0.  +l  -0. )  ]} 


n=0 


F  IN  F 

an  A9  +  -  on  )A0 


F  N 
+  a«+i 


•)  A0 


Po 


v0 


-  S« 

n=0 

^  ,  F  ^  N  > 

4n  A  °r  ar  +x 
n=0  n  ■  **  La 

Afl.  ^  /  F  N  . 

=  "  re.  +Tre.+[ 


(2-nPo) 


„  A0  ^  /  F  N  . 

*  =  417  Zj.  ar  +  ar.+l) 
n=0  »  » 


(2-31Vo) 


(2-31^ 
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<2-31C0) 


In  case  there  are  no  discontinuities,  the  "near"  and  "far" 
values  of  any  loading  stress  are  alike  at  each  division  point  and 
we  can  drop  the  N  and  F  designation  and  obtain  the  following  ex¬ 
pressions  for  the  zero  terms  of  the  loading  coefficients: 
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REPRESENTATION  OF  COEFFICIENTS  FOR  # (z)  and  Y (z) 
IN  TERMS  OF  LOADING  COEFFICIENTS 


Our  next  task  is  to  express  ak  and  a^ ,  the  coefficients  in 
the  series'  for  $ (z)  and  Y  (z)  respectively,  in  terms  of  pk  ,  pk  , 
vk  ,  v*  -  riy  ,  ,  and  ,  which  are  the  Fourier  coefficients 

representing  the  loading  on  the  boundaries.  Since  we  have  found 
these  loading  coefficients  in  the  previous  paragraph  in  terms  of 
the  given  loading  stresses,  our  solution  will  be  formally  com¬ 
plete  when  we  obtain  a,,  and  ak  in  terms  of  pk  ,  ,  vk  ,  etc. 

Again,  because  the  Fortran  system  cannot  allow  negative  sub¬ 
scripts,  we  shall  formulate  all  series'  and  their  coefficients 
to  avoid  such  use.  Therefore,  we  shall  employ  the  following 
equivalence  of  symbols: 

A_k  —  A,,  —  rw  +  iC» 

Alk  =  =  p. k  +  i  v* 

a-k  =  ~  £Lc  + 

a-k  =  *  It  +  i&k 

Thus,  we  must  express  a,  8.  y.  6,  a,  8,  y,  6  in  terms  of  p ,  p,  v, 
T|,  jn,  C# 

Coefficients  of  $ (z) 

From  the  first  of  Eqs.  (2-12)  and  Eq. (h)  we  obtain  a0. 

Oo  =  Re  a„  =  (Po+ivo  (Tb+i-Co)  *  Po  -*i  n>  (2_32a) 

2  (R|  -  R?)  2(R»-R») 

80  =  any  real  number.  Values  for  8o  produce  rigid  body 
rotation  (see  earlier  paragraph  on  displacements) 
and  will  be  assigned  for  convenience  after  a 
specific  problem  is  solved  on  the  computer.  During 
computation.  Bo  *  0.  (2-32a) ' 

We  obtain  ax  and  Bi  from  Eq. (2-20)  after  changing  Alx  to  A/ ,  etc. 
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ax +i3i 


.  *TH|-  A,'»?  .  2A.'R. 

RJ-  R{  (Rj  +  R? )  (x+1) 

=  Rj  (_pt  -ivi )  -R?  (a  -l£i )  .  2  (th  +1Ci)Ri 
RJ  -  R{  {R|+I^ )  (m+1) 

•  „  m  **1  Pi "  R?  %  2rii 

R*-  R*  "  (Rj+R?)  (x+1) 

(2-32b) 

g  *  _  2Ci 

R*-  R{  (R|+R?)(k+1) 

The  remainder  of  the  ak  and  Bk  for  k  ■  2,  3,  . . .  «  are  computed 
from  Eq. (2-21)  with  Dk  given  by  Eq. (2-16) . 

„  _  (1+k)  (R3-  R?)  (A^RJk+a-  A*'  Rf  k  +s  MR^a  k  +a  -  R{*8k+a)  (^Rj+a.  ^'R^8) 

“k  1  1  p 

Dk 

«  lll*Li£b  HI! (Pk^Vk)^'*'3  -  (tw-t-iCk)^4,8] 

Dk 

.  (R^3k  ~*~a -  Rrak+a)CP>-lv>)«S+a  - 

Dk 


•  «(Xk 


0k 


.  (1+kXR3-  R?  Xpk  R^k +a-  r\t  Rfk  +a  Hyt+a’  «Ta  » +a  X^  R$ +a  -  jp,,  R{ +a ) 

°k  (2-32c) 

*  (^XRg-  R? Xv„ R^k +a  -  CkRTk+a)f(Riak+a-  RTak4,aXv)tR5't>a-  ^R?*3) 

Dk 


We  can  determine  ^  and  ^  from  Eq.  (2-18)  and  these  terms  are  given 
below  in  Eqs. (2-33a) 


a-i 


o.+iei  =  *l*l  =  (Jb>z±kl*L 

-*  ~  H  +  l  H+l 


a i 


HiiL 

H  +  l 


& 


-CiRg 

K+l 


(2-33a) 
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The  remainder  of  the  ak and  fik  for  k  *  2,  3,  . . .  ®  are  also  ob¬ 
tained  from  Eqs. (2-21)  and  (2-16)  but  with  k  replaced  by  -k  and 
utilizing  the  changed  symbols.  Note  that  D.k  *  . 


a.k  =  dv  =  Uy+ifiy 

(1-k)  (R|-  R?)  (A,;'R*+2-  A^R**8) 


Sk 


(R%*+8-  r2  k+8  )  (A^  R^*  +a  -  A^Rf**2) 

Dk 

=  (l-kXR|-  Rf  x_pk  R*  +a  -  R?»+8XPyRr+a-  tt.RT^3) 


£k  = 


(2-33b) 

(1-kXRg-  R?  Xvk  R| +2  -  CkR5+aMR?k+a-  Rj k  +a  Xvk  R^k  +a  -  CkRl‘k+a) 


For  efficiency  of  computation,  notice  the  common  parts  in 
Eqs.  (2  32c)  and  (2-33b) . 


Coefficients  of  Y  (z) 

We  may  obtain  yk  and  6k  for  k  =  0,  1,  2,  ...  ®  from  Eq. (2-22b) 
or  Eq. (2-22c) .  As  explained  earlier,  the  use  of  Eq. (2-22b)  re¬ 
sulted  in  loss  of  precision;  therefore  we  shall  use  Eq. (2-22c)  . 

In  some  check  problems,  it  was  verified  that  it  makes  no  difference 
which  one  of  Eq.  (2-22b)  or  Eq.  (2-22c)  is  used  as  long  as  all  com¬ 
ponents  are  accurately  computed. 

=  yk  +i6k  =  -(l-»k)ak+a  +  j^+j-  A£+a  for  k  *  0,  1,  2,  ...  ® 
=  —  (1+k)  (au  +a  +  4-8  )  t  (oty  +a  —  ij3y  +a  )  ~  (Py  +a  +  *-Vk  4.3  ) 


Yy  =  -  (l+k)ak+a+  a*  +a”  Oy+a 


6y  —  “  (1+k)  0k  +8  -  Jjt +a  “  Vy  +a 

We  obtain  Yi  and  from  Eq. (2-19) . 

nAl,R1 


a-i  *  *1  ®  Yi+  =  - 


1+K 


for  k  *  0,  1,  2, 


Mrn+iCiJRi 
- ThJ - 


(2-34) 
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KTfrRt 

1+K 


(2-35a) 


(2-35b) 


The  remainder  of  the  ^  and  6*  are  found  from  Eq.  (2- 22b)  with  k 
replaced  by  -k  and  using  the  change  of  symbols  to  eliminate  nega- 
give  subscripts. 

a-k  *  *£  *  lm+  il*  *  (k-1) Rj  a8 _k  +  ak -aR?k "a-  Ag'_k 
=  (k-l)Rjall_a  +  S-aR?1"8-  A„'-sRi 

*  (k-l)R?  (ak-a  +  i_Bk-a)  +  (ak_a-  i0t -a  > (]> -a+  i£k-a>*i 


•  *.  X*  *  (^-l)R?Sk-a+  a-k  _a  H? 11  -  r^-aRi 

Ak  *  (k-l)R?J^-a-  Bk-aR?k"3-  ik-a^i 


(2-35c) 

for  k  *  3 ,  4 ,  ...  ® 


Again  it  was  found  in  some  check  problems  that  ^  and  6*  could  be 
found  from  Eq.  (2-35c)  based  on  Eq.  (2-22b)  or  they  could  be  found 
from  a  relation  corresponding  to  Eq.  (2-35c)  based  on  Eq.  (2-22c) ; 
both  methods  gave  the  same  values  for  ^  and  6*  if  the  components 
were  accurate.  Here,  however,  it  is  better  to  compute  with  Rt 
included  because  two  terms  in  each  of  Eq.  (2-35c)  contain  1^  to 
positive  exponents  in  k  and  these  two  terms  become  negligible  com¬ 
pared  to  the  other  term  for  large  k. 
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All  that  remains  now  is  to  express  our  stresses  and  dis¬ 
placements  within  the  ring  in  terms  of  series'  using  no  negative 
subscripts. 
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FINAL  EXPRESSIONS  FOR  STRESSES  AND  DISPLACEMENTS 
USING  NO  NEGATIVE  SUBSCRIPTS 


Stresses 


From  Eqs. (2-26) ,  (2-27)  ,  and  (2-28)  we  obtain  the  three 
stress  components  by  splitting  each  series,  substituting  -k  for  k 
in  the  part  summed  from  -®  to  -1,  and  by  employing  the  modified 
symbols  used  above  to  avoid  negative  subscripts. 


o0=  y^rk  |(2+k)  (gk  Cos  k0-04Sin  k0)+Y*Cos(k0+20)-6kSin(k0+2e) 


+  ^  r"k[(2-k)  (a*  Cos  ke+^Sin  k0)+^  Cos  (2q  — ke )- 6^  Sin  (20-k0)  J 

(2-36a) 

GO 

ar  *  ^rk  [(2-k)  (av  Cos  k0-0kSin  k0 )  -yk  Cos  (k0+20)  +6k  Sin  (k0+2g )  J 
00 

+  )  r“k[(2+k)  (a*  Cos  kO+j^Sin  k0)-^  Cos  (20-k0)+j&*Sin(20-k0)  I 
k-1  L  J 

(2-36b) 

00 

t  _  —  £  rk[  k0k Cos  k0+kakSin  k0+6k Cos  (k0+20) +yk  Sin  (k0+20 )  J 


r0 


k=0 


+  ^r”k  j^-k£k  Cos  k9+kakSin  kO+^Cos  (20-k0)+j^kSin(20-k0)  J 

(2-36c) 


Displacements 


We  find  displacements  from  Eqs.  (2-4a)  and  (2-4b)  with  no 
negative  subscripts  by  the  same  sort  of  substitutions  as  were 
done  for  stresses. 

r-y  rl  r  ") 

2pu  *  Coa  *8+1*  sin  *8 )  +[(*£1  -li ) COS  0  +  (kJ3x+6i  )Sin  0Jln  r 

£*ri  +k  r-,4. 

+  Co*  *8-8* Sin  k0)-)rx+k  (akCos  k0-0kSin  k0) 

0  0 
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-  £rl_k  (a*  Cos  kO+g^Sin  k0)  -^[3^[ljtCo8(20-k0)-Ak  Sin(20-k0)j 

2 


u 

1 


-  [yk  Cos  (k@+20 ) -6k  Sin  (k0+20)  J 

0 

2\jl\i  =  [(Haj-^i)Cos  0+(Hgj +6^  )  Sin  ejln  r  -  Cos  0  -  01Sin  0 
00 

+  £rl“k  (*y\~)  (OjcCos  kO+j^Sin  k0) 


or 


00 

+  £r1+k  (”"^-)  (ak Cos  k0-Bk  Sin  k0) 

0 

00  1  k 

-  [^JtCos(k0-20)+6J{Sin(k0-20)  J 


00  ^ 

-  ^YTk  [Yk  Cos  (k0+20 )  -6k  Sin  (k0+20 )  J 


(2-37a) 


2|JV  (”g»sin  k0+l*Cos  k0)  +  [h  (ix  +6j  )  Cos  0  +  (^-HOj)Sin  o]ln  r 


+  (o-n  Sin  k0+0k  Cos  k0)+£rx+k  (akSin  k0+3kCos  k0) 

0  0 

00  00  .  — ^ 

+  ^rx”k  (-a* Sin  k0+^jj  Cos  k0)+^~£^kSin(20-k0)+6J,Cos(20-k0)  J 

®  ^  ^Kk 

+  £ykSin(k0+20)+6kCos(k0+20)  J 
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or 


2|iv  *  £  (k^  +6j,  )  Cos  9  +  (^i  -ho* )  Sin  ejln  r  +  JjCos  0  -  ajSin  0 
00 

+  ^rl“k  (-a* Sin  kO+g^Cos  k@) 

m 

+  ^rl+k  (K-^^k)  (ak  Sin  k0+0kCos  k0) 

0 

®  1  -k 

+  Yt*  [“Ijt  Sin(k0-20)+6Jt  Cos  (k0-20)  ] 

2 

w1  "*"k  r  "i 

+  £-[-£  |_Yk  Sin  (k@+20 )  +6k  Cos  (k0+20 )  J  (2-37b) 

0 

In  Eqs. (2-36a,b, c)  and  (2-37a,b)  we  obtain  a  and  0  from 
Eqs.  (32a, b,c)  ;  a.  and  £  from  Eqs.  (2-33a,b)  ;  y  and  6  from  Eqs. 
(2-34)  ;  and  x  and  6_  from  Eqs.  (2-35a,b,c)  . 


Alternate  Forms  for  Stresses 


While  Eqs. (2-36)  are  quite  compact  in  appearance,  there 
are  some  possible  objections  to  the  form  of  the  expressions: 

(a)  the  trigonometric  functions  of  0  have  three  separate  values 
of  argument,  and  (b)  the  coefficients  y »  X'  JL  mu®t  be 

computed  and  stored. 

We  would  like  to  have  only  one  value  (namely,  k0)  of  the 

trigonometric  argument  in  the  kth  term  so  that  we  have  a  true 
truncated  Fourier  series.  Such  a  series  forms  the  most  accurate 
trigonometric  series  possible  for  k  terms.  It  would  be  possible 
to  make  large  errors  if  only  a  part  of  the  Fourier  coefficient 
of  the  last  term  were  included;  this  could  happen  if  a  large 

missing  part  nearly  canceled  the  included  part  to  make  a  small 
total  coefficient. 

The  reason  we  would  like  to  eliminate  y,  x*  and  6_  is  to 
reduce  storage  and  to  avoid  some  possible  lack  of  precision  in 
their  computation.  If  all  coefficients  are  to  remain  stored  (as 
was  done  in  first  Fortran  version) ,  this  elimination  of  necessary 
storage  may  be  important. 
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Therefore,  in  Eqs.  (2-26),  (2-27)  ,  and  (2-28)  ,  we  shall  re¬ 
assemble  the  various  terms  so  as  to  have  our  series'  expressed 
with  coefficients  of  sin  k6  and  cos  k0;  in  other  words,  the 
stresses  themselves  are  just  Fourier  series'  even  though  the 
Fourier  coefficients  are  involved.  We  shall  also  express  y, 

6,  and  t_  in  terms  of  a,  a,  0,  p,  v#  and  ^  from 

Eqs.  (2-34)  and  two  relations  giving  ^  and  from  Eq.  (2-22c)  . 
These  two  relations  (for  R?  =  1)  are  as  follows: 


Xjc  =  (k-l)a*-2  +  ak  _3  -  £*-3 

Ak  ~  (k-l)^  _2  -  Bit  “3  ”  "3 


00 


(q) 


Omitting  considerable  algebra,  we  obtain  the  stresses  from 
Eqs. (2-36a) ,  (2-36b) ,  and  (2-36c)  in  the  above  indicated 
modified  form  as  shown  by  Eqs. (2-38a) ,  (2-38b) ,  and  (2-38c) . 


0  ■  2a0  +  5~  (.psrH? ) 

8  1-R» 


+  ^(3r+r“*)a1  +  (r-l+2r_3)ai-  - — ^  -  r**g,  Jcos  9 


+  [-(3r+ +  2r~«)g1  +  —  ■K.Rl.6?  -  r-*Vl]sin9 

L  1+H  ”  J 


+  £  {rk[(2  +  k-  kr"8  + 
2 


r-a  +  r-ak-a)ak_ 


+  r~k  (2  -  k+ kr"8+  r8k  ”8+r"8  )ak -r”8gk  J  J-Cos  k0 


+  |~  (-2  -  k+  kr~8-  r“8-  r"8k  "8 )  0k  +  r“8vkJ 

+  r~k  £  (2  -  k+  kr“a+r”8  +  r8k”8)0k-  r^v^Jj-Sin  k0  (2-38a) 


I 
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o  =  2a0  -  r~3R?  (Po-lo  ) 
r  1-Rj 

+  f(r  -  r-»)a,+(3r-1-2r-*)al+  —  xR>  H  +  r-»  Pl  ]cos  e 

1+K  “ 

+  [-(r  -  r-»)B,  +  (3r->-  2r-»)6l-  —-*3  +  r^v.l  Sin  0 

1+H  “  J 

CO 

+  k^2^^<2  "  k  +  kr_3’  r"3'  r'ai,"S)a‘+ 

+  r-‘[(2  +  k  -  r-»-  kr-»-  r>‘-»)a,+  r->£<]  }cos  ke 

ao 

+  Z  {rk[^“2  +  k  -  kr“3+  r"8+  r“2  k  “2  )  9,  -  r~aVkl 

k  —  2  J 

+  1  k[(2  +  k  “  r”3“  kr“a-  r2k”-)gk  +  r”avkJJ-Sin  ke  (2-38b) 


T  =  r“aR?(v0-Co) 


r9 


1-R? 


+  f  (r  -  r-^R}  )Bl  +  (-r-1+  2r-aRa)e1-  — *R*  C*  -  r“«^Ci  J 


1+H 


Cos  9 


+  [(r  -  r^RjJax+tr-1-  2r"«  R?  )  ax  -  £'**%  ^  -  r"3  R?  lsi 

1+H  “  J 

OO 

+  y  {r“  [Bi<  <k  -  *r*a  +  r-a-  r'a*'1-  r*sv,  J 

+  r-‘[s,  (-  k  +  kr-a+  r-a-  r*'-9)  -  r-av,]  }cos  ke 


Sin  e 
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w 

+  £  jrw £ak  (k  -  kr"a  +  r“a-  r"3k  “a )  -  r~apkJ 

2 


+  r”k  £ak  (k  -  kr“a-  r~2+  rak“2)  +  r“a£k  J  |sin  ke  (2-38c) 


In  equations  (2-38a,  b,  c)  the  value  of  Rg  =  1  has  already  been 
used;  similarly,  computations  of  a#  B ,  &,  and  g  should  be  done 
with  Rg  set  equal  to  1. 

For  efficiency  in  computation,  notice  that  in  Eq. (2-38a)  the 
coefficient  of  ak  is  the  negative  of  the  coefficient  0k ;  also  that 
the  coefficient  of  gk  is  the  same  as  that  of  8k .  These  same  re¬ 
marks  hold  for  Eq.  (2-38b)  ,  and  negatively  for  Eq.  (2-38c)  .  In  fact, 
if  the  coefficient  of  ak  in  Eq.  (2-38a)  is  named  XA,  then  the  co¬ 
efficient  of  ak  in  Eq. (2-38b)  is  (4-XA) ;  similarly  if  we  label  the 
coefficient  of  ak  in  Eq. (2-38a)  as  XB,  then  the  coefficient  of  ak 
in  Eq. (2- 38b)  is  (4-XB) . 


Alternate  Forms  for  Displacements 

For  the  same  reasons  as  for  stresses,  we  would  like  to  express 
displacements  with  only  one  argument  in  0  in  the  general  term  and 
also  eliminate  the  use  of  y,  6,  y,  and  6.  Again  we  do  this  by  re¬ 
grouping  the  terms  of  the  several  series  and  substituting  for  y ,  6 , 
6,  their  values  from  Eqs. (2-34) ,  (2-35a,b) ,  and  (q)  .  From 
EqsT(2-37a)  and  (2-37b)  ,  we  obtain  Eqs.  (2-39a)  and  (2-39b)  . 


2(ju  =  r(x-l)a0  + 

1-R? 

+  [(xat  +  ) In  r  +  ra  +  '^2_(2al+  a,-  Oi ) ^jcos  9 

+  [(xB,-  r  -  Bi-  ^-(26!  -  B!  -  v,  )-ra  (*=^)8,  ]sin  9 

» 

+  [r(*-i-k)+  r-1  <i+k)+  r"ak"i]ak 

2 
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+  jrg  [r  (n-l+k)  +  (1-k)  r~l  +  r8k~l  ]ak  -  — £k  +  ^f-pk  }Coa  *9 

00 

+  Z  {■  kTT  [*■  (K-l-V)  +  (l+Tc)  r— 1  +  r-a*-‘]B, 

+  [r (h-1+1c)  +  (1-k) r"1  +  r^~‘ 

2uv  «  r  (k+1)  Bo  ~  —  1.^9_t°)Rl 

1-R? 

+  [(kBi-  r  +  Bi+  ra  (^)  Bi  -  £^-(2B,  -  81-  VjlJcos  6 

+  [-(Kas+  K-^.-R? )  In  r  -  2l+  ra  (^>1,+  ^(25,+  at-  g,)]sin  9 

00  k 

+  I  {f^k  |"r(*+1+k>  +  r-^-Ml+lOr-1  ]ek 
2 

+  jrg  [r(n+l-k)-  r8k  +  (1-k)  r~l  Jpk  +  — vk-  j^-vk}cos  ke 

+  £  {■jjj  [r(n+l+k)  +  r~8  k  “l  -  ( 1+k)  r "l  Jak 
2 

+  jrg  [-r(K+l-k)+  r8k  ~x  -  (1-k)  r”1  Jak  - 


rk  t 

■£k  -  l[riPk}sin  ke 

(2-39b) 


rk-l  <, 

-  k^1"-vk  Jsin  ke 
(2-39a) 


Modified  Expressions  for  Stresses  and  Displacement a  to  Avoid  Over¬ 
flow  of  Computer 

Because  the  separate  computation  of  such  terms  as  r~k  (a  number 
less  them  one  to  a  large  negative  power)  resulted  in  overflows  in 
the  computer,  it  was  necessary  to  modify  Eqs.  (2-38a,b,c)  and 
(2-39a,b)  considerably.  Basically,  the  indicated  expressions  were 
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expanded  so  as  to  perform  algebraically  the  subtraction  of  large, 
approximately  equal  numbers  rather  than  attempt  to  do  it  numerically 
by  computer.  Then  the  various  terms  in  each  expression  were  grouped 
in  such  a  mam.er  that  no  single  term  "ran  away".  This  required  a 
considerable  amount  of  manipulation  which  will  not  be  repeated 
here — only  the  final  expressions  are  shown  below.  Actually,  many 
individual  terms  in  the  expressions  below  tend  to  zero  as  k  becomes 
large,  and  an  underflow  condition  could  result.  Underflows  were 
avoided,  howe/er,  by  programming  tests  prior  to  calculation — then 
if  the  number  to  be  computed  would  underflow,  the  number  was  merely 
set  to  zero  and  not  computed. 

In  the  following  modified  forms  for  Eqs. (2-38a,b,c)  and 
(2-39a,b) ,  we  shall  use  the  following  terminology  to  avoid  excessive 
writings 

Tl  ■  1  +  k  -  Rj  -  kR? 

T2  =  1  -  k  -  R?  +kR* 

H  =  ka  (2R*  -  R*  -  1)  +  R?  (R?k  -  2) 

T9  “  1  +  HR?‘~a 

0L„  =  {Rf-Sl[pk  <T1)-Pk  ]-  R{[ri,(Tl)+  3_Rj]  +  2)IIR?k+  pk}(T9) 

6*  »  {R?‘”a[',k  (Tl)+v„]+  R$[-Ck<Tl)+  C„ ]-  yk}(T9) 

T3  =  2  +  k  -  kr”a  +  r”a 

T4  =  2  -  k  +  kr”a  +  r_a  +  ra  k  ~a 

T5  =  ak  r”k 

=  (|s-)',{R?_a[(T2)ok-  0k  ]-R?k  ((T2)  ry  +  R?  ri,  ]  +  3,  +  pkR?‘}»T9) 

T6  =  r”k  “2  (ak  -£k ) 

=  ^‘^{Rf-'jtTDp,-  Pk-  Hpk]  -  RTJrv(Tl)  -  3,  +  3,R?‘-a}(T9) 
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T7  *  r"k  Bk 

=  (5i-)‘{R{-a(tT2)v,-  vk]  +  Rf*[c,Bf-  Ck  (T2)  ]  -  Ck-  vkRj‘}(T9) 
T8  =  r"»  -»  (-0,  -  vk  ) 

=  ll1-* "  +S  {*»  "*  [vk  -  vk(Tl)-  Hvk]  -  RTaCk(Tl)+  Ck-  CkR?l,',}(T9) 


T10=  k  -  kr”a  +  r”a 

Tll=  -k  +  kr“a  +  r”a  -  rak~a 


T12= 

k+1 


+  r"1 


Tl3= 


r 


k-1 


r-i 


+ 


rak 

k-1 


T14=  r 


U+1+*) 

1+k 


r-i 


From  Eq. (2-38a) , 


2af 


.-a  Da 


(Po- 


1-R? 


+  [  (3r  +  r“*  )  ax  +  (r“x  +  2r“*  )  ttl  -  - — oi  "|co»  e 
+  r.(3r  +  r“®  )  +  (r“l  +  2r"*  )  +  ■-  1  kRi  -  r“®  Vl  lsin  0 

L  “  1+K  ”  J 

ao 

+  £  {rk[(T3)ak-  okr"a]  +  (T4)  (T5)  +  T6  }cos  k0 

1C*  2 

00 

+  £  {r*[-(T3)0k+  r“avk]  +  (T4)  (T7)  +  T8}sin  k9  (2-38aa) 
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From  Eq.  (2-38b) , 


ar  *  2<xc 


_  r~aR?  (p0-  r\p ) 
1-R? 


+  r  (r  -  r^Jat+Or-1-  2r-®)ai+  r  +  r“*  Pl  Icos 

L  “  1+K  “  J 


f-(r  -  r-3  )  Si  +  (3r-1  -  2r-3)g1-  C.liS-k.  +  r^lsi 
L  “  1+K  -  J 


Sin  0 


+  £  {rk[(4-T3)ak  +  r”apk]  +  (4-T4)  (T5)  -  T6}cos  k0 
k=  2 

00 

+  £  {rk[(T3-4) 0k-  r"avk]  +  (4-T4)  (T7)  -  T8|Sin  kQ  (2-38bb) 


From  Eq. (2- 38c) , 


Tr9  .  r-»R?(v0-C0) 

1-B? 

+  [<r  -  Rj  )  S,  +  (-r“l  +  2r-«R?)Bl-  r"1  '  r**I?£i]co8  8 

+  [(r  -  r-»R{  Jm  +  tr-1-  2r-‘-R?)al  -  — KR>  "i-  -  r"3  R?  n,  ]sin  9 

L  "  1+K 

r"3 vk ]  +  (Til)  (T7)  +  T8}cos  k0 
r“a  pk ]  -  (Til)  (T5)  -  T6}sin  k0  (2-38cc) 


+  Y  {rk[ (TlO) 0k  - 
k«2 L  L 


+  y  {rk[(TlO)ak- 
1c®  2 
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Prom  Eq. (2-39a)  , 


2UU  =  r(K-l)a0  +  flliS  iSa^Jai 

l-R? 

+  ^2i+  "  1  ln  r  "  5‘  +  £T'(22»+  «i-  Oi)  +  ra  (^)ai]cos  8 

+  ^K-,_  * ln  r  "  "  •£2~<2§x  ~  6i-  vj  -  r»  (*=£)  B,  ]sin  0 

+  £  {rk  (T12)ak-(T13)  (T5)  +  ^ +  £i2|I}cos  k8 

+  £  {-*k  (Tl2)0k-(T13)  (T7)  -  — 1  v*  +  kQ 

k=2 L  k-1  k+1  J 

(2-39aa) 

From  Eq. (2-39b) , 

2 pv  *  r(H+l)0o  -  — — ^V°  ~ 

l-R? 

+  C**-1"  ~u!r>ln  r  +  2»  +  rSt^)Bi  -  0i-  vjjJcos  e 

+  [-(KOX+  K^%1n  r  -  2l+  +  ^(2^  +  ax-  gi )  ]sin  8 

+J2{r‘  (T14)e*  +  <T13  -  £T>  <T7>-  '  i^T<T8)}cos  ke 

00 

+  y  {**  (T14)a,  -  (T13  -  j^)  (T5)  -  fill  p»  +  j^j<T6)}$in  ke 
k=2  k- X 

(2-39bb) 
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CHAPTER  2 


THE  COMPUTER  PROGRAM  AND  ITS  USE 


General 

The  FORTRAN  II  computer  program  for  solution  of  the  thick  cy¬ 
lindrical  shell  by  the  Muskhelisvili  method  will  be  described  in 
such  a  manner  that  it  can  be  used  by  persons  quite  unfamiliar  with 
the  theoretical  development. 

Basically,  the  program  will  yield  stresses  and  displacements 
at  all  specified  locations  in  the  shell  for  any  specific  stress 
loading  on  the  boundaries.  This  loading  must  place  the  shell  in 
static  equilibrium.  Only  plane  problems  can  be  solved;  the  load¬ 
ing  stresses  must  not  vary  in  the  direction  of  the  cylinder  axis. 
Either  or  both  boundaries  may  be  loaded  by  radial  stresses,  shear 
stresses,  or  by  both  radial  and  shear  stresses  (in  any  case,  over¬ 
all  static  equilibrium  must  exist) .  There  need  not  be  symmetry 
in  any  of  the  loading  stresses.  The  loading  stresses  may  be  con¬ 
tinuous  or  discontinuous. 

The  analysis  is  based  on  approximation  of  the  distributed 
loading  stresses  by  straight  line  segments  between  equally  spaced 
points  (called  division  points) .  Therefore,  to  prepare  a  problem 
to  be  solved  by  this  program,  the  ring  is  first  divided  angularly 
into  some  number  (M)  of  equal  parts.  Next  the  numerical  value  of 
each  loading  stress  at  each  division  point  is  divided  by  the  larg¬ 
est  absolute  value  of  any  stress  on  any  boundary.  These  stress 
ratios  at  each  division  point  are  the  loading  stress  inputs  to 
the  program.  Also,  the  program  takes  the  outer  radius  to  be  unity 
and  any  other  radius  as  the  geometrically  similar  fraction  thereof. 
Thus,  there  is  a  bit  of  dimensional  analysis  in  first  reducing  the 
actual  problem  to  proper  computational  form  and  then  interpreting 
the  results  in  terms  of  the  actual  problem.  In  all  that  follows 
here  we  are  referring  to  these  "reduced"  or  "program"  stresses 
and  dimensions. 
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THE  PROGRAM  PROPER 


Brief  Functional  Outline 


Initially  the  problem  and  control  parameters  are  read  in 
then  the  loading  stresses  on  both  boundaries  ("far"  points)  are 
entered.  If  any  discontinuities  exist,  the  "near"  values  of  the 
stresses  at  all  discontinuity  points  are  read  in.  Next  the  loca¬ 
tions  at  which  stresses  said  displacements  ?re  to  be  computed  are 
read.  Certain  boundary  points  are  designated  test  poir  s  at  which 
locations  the  computed  radial  and  shear  stresses  must  match  the 
corresponding  loading  stresses  within  a  certain  tolerance — this 
is  the  fundamental  criterion  upon  which  is  based  the  decision  to 
proceed  with  more  computation  or  to  consider  the  problem  solved 
and  print  results.  The  final  data  input  is  a  table  of  sines  and 
cosines  of  the  angle  for  each  division  point. 

In  case  only  one  boundary  is  loaded,  certain  loading  coef¬ 
ficients  are  zero.  When  such  a  loading  condition  exists,  those 
loading  coefficients  are  set  to  zero  to  avoid  wasted  computation 
time. 


Since  the  entire  method  is  based  on  static  equilibrium  of 
the  shell,  checks  are  made  for  this  condition.  If  the  shell  is 
not  in  equilibrium,  the  program  pauses  and  prints  such  a  message. 
These  static  checks  are  made  after  the  loading  coefficients  for 
k  =  0  (moment  check)  and  k  =  1  (force  check)  are  computed.  Thus, 
very  little  time  is  wasted  if  static  equilibrium  does  not  exist. 

As  soon  as  a  set  of  loading  coefficients  is  computed  for  any 
k  2  1,  complete  use  is  made  of  this  set  of  coefficients  and  they 
are  never  used  again.  Then  the  next  set  of  loading  coefficients 
are  stored  in  the  same  locations.  The  contributions  to  stresses 
and  displacements  by  each  set  of  loading  coefficients  are  accumu¬ 
lated,  so  we  have  a  running  set  of  results  at  all  times.  This 
scheme  has  both  advantages  and  disadvantages.  The  disadvantage 
is  that  stresses  and  displacements  can  be  found  only  for  that  set 
of  computation  points  originally  read  in.  The  advantage  is  that, 
since  only  eight  storage  locations  are  used  for  the  loading  co¬ 
efficients  (for  k  *  1) ,  there  is  no  limit  as  to  the  number  of 
terms  (value  of  k)  that  can  be  used  to  solve  a  problem  having 
complicated  loading  stress  distributions.  Computer  storage  will, 
of  course,  limit  the  number  of  computation  points;  but  this 
usually  will  not  be  a  problem.  If  results  at  more  points  are 
desired,  run  the  program  again  for  the  new  points. 
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The  Flow  Chart 


6 
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O] 


Read  Locations  where  Stresses  and  Dis¬ 
placements  to  be  Computed 


Read  R(JJ),  TH(JJ) 


Read  Computation  Point  No.,  Boundary  No., 
and  Division  Point  No.,  at  which  Loading 
Stress  tests  to  be  made 
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99 


IF(K-l)/  + 


Q 
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0 


IF (K-KINT) 


TOL*TOL+DELTOL 


K-K+l 


^  Is  \ 
K  up  to  KMAX? 
IF (K-KMAX) 

OfYes 


KTES=»-1 


error 


PAUSE 


Print:  Problem 
Complete 


/  Does 
lother  Problem 
Follow?  ^ 


no 


DIMENSION  XLS(?»102)tXLT(2» 102).SN(102> ,CS( 102* , JDI SC U ) , XLSN( 2*4) 
1»XLTN(2»4),NT< 12),R( 12»,TH( 12) ,SIGR( 12) .SIGH  12) » T AU< 12 ) , U ( 12 ) , V ( 1 
22 ) » LL ( 12 )  *J  l  (  12) 
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INPUT  PARAMETERS  AND  DATA 


Definition  of  the  Numbers  in  the  Dimension  Statement 


XLS (A , B) ,  XLT (A, B)  :  A 

B 

SN  (A)  ,  CS  (A)  :  A 

JDISC (A) :  A 

XLSN (A, B) ,  XLTN (A, B) :  A 

B 

NT  (A),  LL (A) ,  JI (A) :  A 

R (A) ,  TH (A) ,  SIGR(A) , 
SIGT (A) ,  TAU (A) ,  U  (A)  , 

V  (A)  :  A 


2  (two  boundaries) 

M  +  2 

M  +  2 

NDISC  (number  of  discontinuities) 

2  (two  boundaries) 

Number  of  discontinuities 

NTEST  (number  of  test  points) 

NRTH  (number  of  computation  points) 


Definition  of  Problem  and  Control  Parameters 


M  — number  of  equal  sectors  into  which  ring  is  divided,  must 
be  divisible  by  4. 

TOL  — tolerance  within  which  all  computed  stresses  must  match 
loading  stresses  at  test  points. 

DELTOL — increase  in  tolerance  after  KINT  terms  of  Fourier  series 
representing  loading  stresses  have  been  used. 

KINT  --number  of  terms  at  which  TOL  will  be  increased  by  DELTOL. 

KMAX  — maximum  number  of  terms  the  program  will  be  allowed  to  run. 

NRTH  —the  number  of  computation  points;  that  is,  the  total  num¬ 
ber  of  locations  at  which  stresses  and  displacements  are 
to  be  computed  including  the  test  points  on  boundaries. 

NDISC  — number  of  division  points  at  which  there  is  a  loading 
stress  discontinuity  of  any  kind  on  either  boundary. 

Even  if  only  one  stress  on  one  boundary  has  a  discon¬ 
tinuity  at  a  given  division  point,  it  counts  as  one  of 
NDISC. 

NTEST  — number  of  test  points.  Test  points  are  those  locations 
on  the  boundaries  at  which  loading  and  computed  stresses 
are  compared.  The  test  points  are  included  in  the  com¬ 
putation  points. 
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LC  — if  LC  ■  0,  we  do  not  print  loading  coefficients; 
if  LC  -  1,  we  do  print  loading  coefficients. 

IPR  — if  IPR  =  0 ,  we  do  not  print  accumulated  stresses  and  dis¬ 
placements;  if  IPR  =1,  we  do  print  them. 

STAT  --amount  within  which  the  computed  checks  for  static  equili¬ 
brium  must  fall.  Since  there  are  various  unavoidable 
inaccuracies,  STAT  must  not  be  zero,  or  the  program  will 
always  stop  because  of  apparent  non- ecru:' librium.  Judge¬ 
ment  in  choosing  STAT  must  be  based  on  total  forces  and 
moments  on  the  ring.  Perhaps  .1%  of  absolute  area  under 
loading  stress  curves  might  be  tried.  If  a  problem  has 
symmetric  loading,  use  a  rather  small  STAT — this  gives  a 
check  on  data  entry  error. 

Rl  — radius  of  inner  boundary  (for  outer  boundary  =  1) .  R1  is 

the  ratio  of  inner  to  outer  radius  in  the  actual  problem. 

RPI  *  1/tt 

DELT  =  2tt/M,  the  angle  in  radians  between  division  points. 

RDELT  =  1/DELT 

XKAP  =  h  =  a  Lame '  constant,  k  =  3-4a  for  plane  stress,  h  = 
for  plane  strain,  where  o  =  Poisson's  ratio. 

LLOC  — if  LLOC  =  -1,  inner  boundary  only  is  loaded; 
if  LLOC  =  0,  both  boundaries  are  loaded; 

if  LLOC  =  +1,  outer  boundary  only  is  loaded. 

UF  — underflow  control  parameter.  If  certain  numbers  become 

less  than  10"^,  we  set  that  number  to  zero  to  avoid  an 
underflow;  of  course,  the  test  is  made  before  computation. 
UF  *  20.  has  been  used  with  IBM  7040  (which  underflows  at 
10"*8),  but  the  characteristics  of  the  machine  being  used 
should  govern. 


Loading  Stresses  and  Their  Specification 

The  loading  stress  notation  is  shown  for  radial  stress  on  the 
outer  boundary.  Similar  notation  is  used  for  the  other  three  load¬ 
ing  stresses.  The  division  point  numbering  starts  with  1  at  9  *  0 
because  a  subscript  cannot  be  zero  in  a  Fortran  subscripted  variable. 
Shown  in  Fig.  16  are  the  radial  stresses  at  several  points  with  no 
discontinuities  and  one  point  at  which  a  discontinuity  exists. 

Only  one  value  for  each  of  the  four  loading  stresses  is  read  in  if 
there  are  no  discontinuities  at  that  division  point  in  any  of  the 
four  stresses.  However,  if  one  or  more  of  the  four  stresses  has  a 
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1 


XI. S  (1, 


i 


Figure  16.  Notation  used  with  loading  stresses. 
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discontinuity  at  a  given  division  point,  all  four  stresses  are  con¬ 
sidered  to  be  discontinuous,  and  two  values  are  read  in  for  each 
stress  (even  though  the  two  values  may  be  identical  for  some  of  the 
stresses) .  Definitions  for  loading  stresses  follow: 


XLS (L,N) 


XLT (L,N) 


N 


--radial  loading  stress  (cr)  on  boundary  L  at  division 
point  N.  If  there  is  a  discontinuity  at  this  division 
point,  XLS (L,N)  is  the  "far"  point  as  shown  in  Fig.  16. 

—shear  boundary  stress  (Trg)  on  boundary  L  at  division 
point  N.  Again,  this  is  "far"  point  if  there  is  a 
discontinuity. 

— division  point  number  starting  with  1  and  proceeding 
through  K  +  2.  This  overlap  of  two  on  the  numbering 
was  necessary  because  of  starting  and  ending  some  of 
the  series' . 


L 


—boundary  identification.  L  -  1  signifies  inner  boundary; 
L  *  2  means  outer  boundary. 


XLSN (L, JJJ) — "near"  value  of  radial  loading  stress  at  a  discontinuity 

on  boundary  L.  JJJ  is  the  serial  number  of  the  discon¬ 
tinuity  starting  with  1  for  the  first  discontinuity 
and  running  through  NDXSC  for  the  last. 


XLTN  (L, JJJ) — "near"  value  of  shear  loading  stress  at  a  discon¬ 
tinuity  on  boundary  L.  JJJ  same  as  above. 

JDISC (JJJ)  — division  point  number  at  which  discontinuity  serial 

number  JJJ  occurs. 


Computation  Point  Identification 

R (JJ)  — radius  of  computation  point  serial  number  JJ.  JJ  starts 
with  1  and  runs  through  NRTH. 

TH ( JJ) — angle  measured  from  0  =  0  of  computation  point  serial 
number  JJ. 

See  Fig.  17  for  sketch  of  computation  point  numbering  system. 


Test  Point  Identification 


NT (I) — computation  point  number  of  the  test  point  having  serial 
number  I.  I  starts  with  1  and  runs  through  NTEST. 


LL(I) — boundary  identification  of  test  point  having  serial  number  I. 

LL (I)  =  1  for  inner  boundary;  LL(I)  -  2  for  outer  boundary. 

jl (i)- -division  point  number  of  test  point  having  serial  number  I. 

Figure  17  shows  the  numbering  scheme  for  confutation  points 
and  test  points.  We  start  with  computation  point  number  1  at  the 
smallest  angle  and  the  smallest  radius  at  which  stress  and  dis¬ 
placement  computations  are  to  be  made.  Then  we  proceed  counter¬ 
clockwise  at  this  same  radius  until  all  computation  points  on  that 
circle  are  numbered.  Next  we  move  to  the  next  larger  radius  at 
which  any  computations  are  to  be  made  and  again  proceed  counter¬ 
clockwise  until  all  computation  points  at  this  radius  are  numbered. 
This  process  continues  until  all  computation  points  are  numbered 
serially  from  1  through  NRTH.  Computation  points  will  be  used  as 
test  points— thus  computation  point  8  may  be  test  point  2.  Test 
points  must  lie  on  division  points  in  order  to  have  a  stress  to 
compare  with.  Test  points  must  not  be  points  of  discontinuity. 

As  an  example,  for  M  =  32,  we  will  locate  on  Fig.  17  the  confuta¬ 
tion  and  test  points  indicated  in  the  following  table: 


Comp.  Pt.  No, 


Test  Pt.  No, 


Radius 

(R(JJ)*fract. 


Angle 


(JJ) 

(I) 

of  Rl) 

(TH  (JJ)  in  : 

1 

— 

.5 

.3 

2 

(1) 

.5 

6tt 

32 

(2) 

16tt 

3 

.  5 

32 

(3) 

42tt 

4 

•  5 

32 

5 

— 

.65 

.3 

.65 

6tt 

6 

32 

7 

— 

.8 

0 

8 

— 

1.0 

.3 

(4) 

4tt 

9 

1.0 

37 

10 

(5) 

1.0 

8tt 

32 
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1,  2,  •••  10  are  computation  point 


numbers 

(1) ,  (2) ,  • • •  (5)  are  test  point 

numbers 

Figure  17.  Example  of  numbering  system  for  computation 

points  and  test  points 

Table  of  Sines  and  Cosines 

In  much  of  the  program,  sines  and  cosines  of  the  division 
point  angles  are  used  many  times.  Since  there  are  only  M+2  of 
these  angles,  the  values  of  sine  =  SN(I)  and  cose  *  CS(I)  are 
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computed  by  a  preliminary  program;  these  values  are  then  read  in  as 
data.  This  procedure  makes  the  program  run  considerably  faster 
than  it  would  if  sines  and  cosines  were  determined  as  functions. 

SN(I)~ sine  of  the  angle  whose  division  point  number  is  I.  I 

starts  at  1  for  6  =  0,  is  M  for  8  «  2w  -  DELT,  is  M+l  for 
0  =  2tt  and  ends  at  M+2  for  0  *»  2tt  -»-  DELT. 

CS  (I)—  cosine  of  the  angle  whose  division  point  number  is  I. 
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OUTPUT  RESULTS 


The  outputs  are  nearly  self  explanatory  in  that  they  will  all 
have  header  symbols  above.  However,  for  clarity  let  us  discuss 
briefly  each  of  the  outputs. 


Stresses  and  Displacements  at  Computation  Points 

These,  of  course,  are  the  results  that  the  entire  program  is 
designed  to  produce.  These  results  may  be  printed  only  once,  after 
the  computed  stresses  match  the  loading  stresses  at  all  test 
points;  if  this  is  desired  set  1PR  =  0.  Another  choice  is  to  print 
these  results  after  the  stress  and  displacement  accumulations  are 
made  for  each  value  of  k;  if  this  is  desired,  set  1PR  =  +1.  The 
latter  mode  enables  one  to  observe  rate  of  convergence,  proximity 
of  computed  to  test  stresses,  etc.  Following  is  an  example  of  the 
form  for  these  results: 


K  R  THETA  SIG  THET  SIGMA  R  TAU  U  V 

38  .65  1.45  1.7  .8  -.9  .1  -.2 


K  =  38 

R  *  .65 
THETA  =  1.45 


•38  Fourier  terms  of  loading  stress  representation 
have  been  computed  and  their  contributions  to 
stresses  and  displacements  accumulated. 

•this  computation  point  has  radius  =  .65. 

•this  computation  point  located  at  angle  1.45  rad. 


SIG  THET  =  1.7 — transverse  stress  (a  )  accumulated  to  date  is  1.7. 

Q 

SIGMA  R  *  .8  ••radial  stress  accumulated  to  date  is  .8. 


TAU  =-0.9  — shear  stress  (t  „)  accumulated  to  date  is  -.9. 

r0 

U  =  .1  —radial  displacement  (u)  accumulated  to  date  is 

u  =  U/2|i  =  .1/2 where  u  =  the  second  Lame' 
constant  (m  =  shear  modulus) . 

V  *  -.2  — transverse  displacement  accumulated  is  v  *  V/2|j  = 

2/2(i . 

Note  that  the  above  example  gives  output  for  only  one  compu- 
tation  point;  if  more  points  were  computed,  the  heading  would  not 
be  repeated  for  each  point. 
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Loading  Coefficients 


It  may  be  desired  to  observe  the  values  of  the  loading  coef¬ 
ficients  ETA, - ZNUB  after  each  k  to  see  how  convergence  is  pro¬ 

gressing.  These  loading  coefficients  are  not  simply  the  Fourier 
coefficients  of  the  loading  stress  functions,  but  are  simple  com¬ 
binations  of  them.  If  it  is  not  desired  to  print  these  coefficients 
after  each  k,  set  LC  *  0;  if  printing  is  desired,  set  LC  *  +1.  This 
output  is  of  the  following  form: 


THE  NEXT  8  NUMBERS  ARE  ETA,  ETAB,  ZETA,  ZETAB,  RHO,  RHOB, 
ZNU,  ZNUB 

0.1  0.2  0.  0. 

-0.3  0.  0.25  -0.2 

Here  we  have  ETA-0 . 1 ,  ETAB=0 . 2 ,  ZETA=0 . ,  ZETAB-0 . ,  RHO*-0 . 3 , 
RHOB=0 . ,  ZNU=0 . 25 ,  and  ZNUB*-0.2. 


Static  Equilibrium  Check 

In  case  the  absolute  value  of  any  of  STCHX,  STCHY,  or  STCHM 
is  greater  than  STAT,  the  following  statement  will  be  printed  and 
punched : 


|  LOADING  STRESSES  NOT  IN  OVERALL  EQUILIBRIUM  | 

After  this  statement  is  punched  and  printed,  the  computer  will 
pause.  Restart  if  another  problem  follows. 

If  a  careful  check  of  the  loading  stresses  reveals  no  error, 
the  above  statement  may  merely  mean  that  the  value  of  STAT  was 
chosen  too  small  for  the  inherent  errors  and  approximations  made. 


Insufficient  Terms  Used 


In  case  the  maximum  desired  number  of  terms  (KMAX)  has  been 
reached  by  k  and  the  computed  stresses  at  the  test  points  have 
still  not  come  within  tolerance  (TOL  +  DELTOL) ,  the  following 
statement  will  be  printed  and  punched: 


MAX  NUMBER  TERMS  COMPUTED  BUT  DID  NOT  MEET  TEST 
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After  the  above  statement  is  punched  and  printed,  the  stresses  and 
displacements  accumulated  to  date  will  be  printed  after  which  will 
be  printed  the  end-of -problem  statement  described  in  the  next  para¬ 
graph  . 


End  of  Problem  Statement 


If  the  problem  is  terminated  because  of  successful  solution 
or  because  of  k  reaching  KMAX  the  following  statement  will  be 
printed: 

|  THIS  PROBLEM  COMPLETE  RESTART  IF  OTHER  PROBLEM  FOLLOWS  | 

The  above  statement  will  not  be  printed  if  the  problem  is  terminated 
because  of  failure  to  meet  static  equilibrium. 

After  the  above  statement  is  printed,  the  machine  will  pause 
and  should  be  restarted  if  another  problem  follows. 


Logic  Errors 

In  case  some  control  number  is  set  to  a  meaningless  value 
accidentally  the  following  statement  will  be  printed: 


|  LOGIC  ERRORS  CHECK  CONTROL  PARAMETERS  | 

If  this  statement  is  printed,  examine  KINT,  KMAX,  LC,  IPR,  LLOC. 
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EXAMPLE  PROBLEM 


Actual  Problem 


A  cylinder  is  loaded  in  plane  strain  conditions  by  compres¬ 
sive  radial  stresses  of  5000  psi  over  two  small  arcs  on  the  outer 
boundary  as  shown  in  the  figure.  The  outside  radius  is  100  inches; 
the  inside  radius  is  50  inches.  Poisson's  ratio  of  the  material 
is  0.3  and  the  shear  modulus  is  2  x  106  psi.  The  loaded  small  arcs 
are  4rr/100  radians.  See  Example  Figure  18. 


Program  Problem  (See  Example  Figure  19) 

Desire  to  compute  answers  at  indicated  locations: 


No. 

r  (inches) 

0  (rad) 

1 

50 

0 

2 

50 

.62831853 

3 

50 

1.5707963 

4 

50 

4.7123890 

5 

60 

.62831853 

6 

62.5 

.62831853 

7 

75 

.62831853 

8 

85 

.62831853 

9 

10C 

0.0 

10 

100 

.62831853 

11 

100 

1.5707963 

12 

100 

4.7123890 

Input  Data 

1.  DIMENSION  XLS( 2,102 ) ,  XLT(2,102),  SN(102),  CS(102),  JDISC (4)  , 
XLSN(2,4),  XLTN (2,4) ,  NT(6),  R(12) ,  TH(12),  SIGR(12) , 

SIGT  (12)  ,  TAU(12)  ,  U ( 12)  ,  V(12)  ,  LL(6)  ,  Jl(6) 
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Example  Figure  18.  Loading  stresses  on  actual  shell 
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1  psi 


1M  (5) 


Poisson's  ratio  =  o  *  .3 


Shear  Modulus  = 
H  =  2  x  10*  psi 


3  (2! 


2\T6 


.62831853  rad, 
\  _  2  Ott 

l  ~  100 

-i — N*=l ,  101 


3i(4J 


(6)  12 


Lutr 

-~l  I  Hfe  r“ 

\ —  N=76 


N=75 


XKAP  =  k 


_  3-o_  =  3- . 3 
l+o  =  1.3 

=  2.0769231 


Let  M  *  100  (each  loaded 
arc  =  2DELT) 

Then  DELT  = 

100 

=  .062831853  rad 
RDELT  =  15.9154943 


Example  Figure  19.  Loading  stresses  on  program  shell 


.2' 


READ  1, 
IPR, 

M,  TOL,  DELTOL 
STAT,  Rl,  RPI, 

,  KINT,  KMAX,  NRTH,  NDISC,  NTEST,  LC, 

DELT,  RDELT,  XKAP,  LLOC 

2  cards: 

Remarks : 

M  - 

100 

shell  divided  into  100  equal  segments 

TOL  - 

.03 

DELTOL  - 

>  .1 

KINT  « 

10 

increase  TOL  by  DELTOL  when  K  «  KINT 

KMAX  » 

100 

end  problem  if  K  *  KMAX 

NRTH  * 

12 

there  are  12  computation  points 

NDISC  * 

4 

there  are  4  discontinuities 

NTEST  - 

6 

there  are  6  test  points 

LC  « 

1 

print  loading  coefficients 

IPR  - 

1 

print  accumulated  stresses  and  displace¬ 
ments  after  each  K 

STAT  * 

.001 

symmetric  problem  should  allow  small  STAT 

Rl  * 

.5 

RPI  = 

.31830990 

reciprocal  of  tt 

DELT  = 

.06283185 

2tt  radians  divided  by  100 

rdelt  * 

15.9154943 

reciprocal  of  DELT 

XKAP  = 

2.0769231 

3— c 

h  =  —  for  plane  strain,  and  c  *  Pois¬ 

son's  ratio  *  .3 

LLOC  = 

+  1 

loading  stresses  on  outside  boundary  only 

READ  1111,  UF 

UF  = 

20. 

READ  2, 
MP2  * 

( (XLS  (L,N)  ,  XLT  (L,N)  ,  N  *  1,  MP2)  ,  L  -  1,  2)  Note: 

M+2  =  102 

204  cards: 

102  blank  cards  for  unloaded  inner  boundary,  L  = 
24  blank  cards  for  L  *  2  and  N  *  1,  2,  ...  24. 

2  cards  XLS(2,N)  XLT(2,N): 


-1.0  0.0  N  *  25 

-1.0  0.0  N  *  26 


1. 
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48  blank  cards  for  L  *  2  and  N  *  27,  28,  ...  74 

2  cards  XLS(2,N)  XLT ( 2 , N ) : 

-1.0  0.0  N  =  75 

-1.0  0.0  N  =  76 


26  blank 

cards  for  L 

»  2  and  N 

-  11,  IB, 

•  •  • 

102. 

READ  3,  (  (JDISC (JJJ) ,  XLSN(L, JJJ) , 

NDISC) ,  L  =  1 ,  2) 

XLTN (L, JJJ),  JJJ  *  1, 

8  cards: 

Re. 

,ks 

• 

• 

JDISC (JJJ) 

XLSN (L, JJJ) 

XLTN (L, JJJ)  L 

JJJ 

25 

0. 

0. 

1 

1 

27 

0. 

0. 

1 

2 

75 

0. 

0. 

1 

3 

77 

0. 

0. 

1 

4 

NDISC  =  4 

25 

0. 

0. 

2 

1 

27 

-1.0 

0. 

2 

2 

75 

0. 

0. 

2 

3 

77 

-1.0 

0. 

2 

4 

NDISC  =  4 

Notes  These  cards  contain  the  " far”  point  loading  stresses  at  the 
four  discontinuities.  See  Example  Figure  19. 


READ  2, 

(R(JJ) ,  TH (JJ) , 

JJ 

=  1, 

NRTH) 

12  cards 

• 

• 

Remarks 

: 

R  (JJ) 

TH  (JJ) 

JJ 

0.5 

0. 

1 

also 

test 

pt.  #1 

0.5 

.62831853 

2 

0.5 

1.5707963 

3 

also 

test 

pt.  #2 

0.5 

4.7123890 

4 

also 

test 

pt.  #3 

0.6 

.62831853 

5 

0.625 

.62831853 

6 

0.75 

.62831853 

7 

0.85 

.62831853 

8 
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1.0 

1.0 

0. 

.62831853 

9 

10 

also 

test  pt.  #4 

1.0 

1.5707963 

11 

also 

test  pt.  #5 

1.0 

4.7123890 

12 

also 

test  pt.  #6 

NRTH  -  12 

See  Example  Fig.  19  for  location  of  these  computation  points. 
These  numbers  are  not  in  parentheses. 


7. 


READ  7, 

(NT (I), 

LL(I),  J1(I), 

1*1,  NTEST) 

6  cards: 

Remarks : 

NT  (I) 

LL(I) 

J1(I) 

I 

1 

1 

1 

(1) 

3 

1 

26 

(2) 

4 

1 

76 

(3) 

9 

2 

1 

(4) 

11 

2 

26 

(5) 

12 

2 

76 

(6)  NTEST  -  6 

See  Example  Fig.  19  for  location  of  these  test  points.  Test 
point  numbers  are  in  parentheses;  e.g.  (4)  . 


Output  Results 

An  example  of  the  output  for  k  ■  1  is  included  here  to  show 
the  actual  format.  Such  intermediate  results  do  not  need  to  be 
printed.  Also  the  printing  of  ETA,  ...  ZNUB  need  not  be  printed. 


1=K  THE  NEXT  8  NUMBERS  ARE  ETA ,  ETAB ,  ZETA ,  ZETAB ,  RHO ,  RHOB ,  ZNU ,  ZNUB 


0. 

0. 

0. 

0. 

0. 

0. 

-0. 

-0. 

K 

R 

THETA  SIG  THET 

SIGMA  R  TAU 

U 

V 

1 

0.50000 

0.  -0.106667 

0.000000  -0. 

-0.  *>41025642 

0. 

1 

0.50000 

0.62832  -0.106667 

0.000000  -0. 

-0.041025742 

0. 

1 

0.50000 

1.57080  -0.106667 

0.000000  -0. 

-0.041025742 

0. 

1 

0.50000 

4.71239  -0.106667 

0.000000  -0. 

-0.041025642 

0. 

1 

0.60000 

0.62832  -0.090370 

-0.016296  -0. 

-0.039452992 

0. 
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1 

1 

1 

1 

1 

1 

1 


0.62500 

0.75000 

0.85000 

1.00000 

1.00000 

1.00000 

1.00000 


0.62832 

0.62832 

0.62832 

0. 

0.62832 

1.57080 

4.71239 


-0.087467 

-0.077037 

-0.071788 

-0.066667 

-0.066667 

-0.066667 

-0.066667 


-0.019200  -0, 
-0.029630  -0, 
-0.034879  -0, 
-0.040000  -0. 
-0.040000  -0, 
-0.040000  -0. 
-0.040000  -0. 


-0.039202052  0. 
-0.039316240  0. 
-0.040096532  0. 
-0.042051283  C. 
-0.042051283  0. 
-0.042051283  0. 
-0.042051283  0. 


Conversion  of  Results  Back  to  Original  Problem 

From  dimensional  analysis  we  find  that  we  can  obtain  a  given 
stress  by  the  rollowing  relation: 


fc} 


max 


=  §-} 


actual  shell 


max 


program  shell 


.  .  S 


(s  ) 

max  act. 

act  ~  (S  ) 

max  prog 


x  S 


prog, 


where 


S  ■  any  stress  in  actual  shell. 

3Ct 

Sp^g  =  corresponding  stress  in  program  shell. 

(S  )  .  =  maximum  absolute  value  of  any  loading  stress  on 

iH3X  3C  u  .  •  •  •  « 

actual  shell. 

(S  )  =  maximum  absolute  value  of  any  loading  stress  of 

max  prog  ... 

program  shell. 


For  our  example. 


(S  ) 
max  act 

(S  ) 
max  prog 


5000 


5000 


Thus  any  stress  in  the  actual  shell  can  be  obtained  by  multiplying 
its  counterpart  in  the  program  shell  by  5000. 


Also,  we  find  from  dimensional  analysis  the  following  rela¬ 
tion  for  displacements: 

fr}  t  -  h  n  =  {* }  . 

^  actual  shell  *  program  shell 


*9 


act 


act 


x  d 


prog 


prog 
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where 


act 

d 

prog 
^  act 


either  displacement,  u  or  v,  in  actual  shell. 

»  corresponding  displacement  in  program  shell. 

:  radius  of  outer  boundary  of  actual  shell. 

Rj pr0g  ■  1  ■  radius  of  the  outer  boundary  of  program  shell. 

For  our  example  problem,  Ra  .  ■  100".  Thus,  any  displacement 

act 

the  actual  shell  may  be  found  by  multiplying  the  corresponding 
displacement  in  the  program  shell  by  100. 


Of  course,  the  usual  similitude  requirements  must  be  met 
before  the  above  conversion  relations  are  valid.  The  shells  must 
be  geometrically  similar,  the  points  under  consideration*1  must  be 
located  geometrically  similarly,  the  elastic  constants  must  be  the 
same,  and  loading  stresses  must  be  similar  throughout.  All  these 
requirements  are  met  by  our  problem  procedure. 
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SECTION  III 


STRESSES  IN  THICK-WALLED  CYLINDERS  WITH  CAPPED  ENDS 

CHAPTER  1 
INTRODUCTION 


This  section  is  concerned  with  static  stresses  in  thick- 
walled  cylinders  with  end  caps.  The  cylinders  are  assumed  to  be 
loaded  axisymmetrically;  thus  the  general  problem  is  to  solve 
axisymmetric  elasticity  problems  involving  a  difficult  geometry. 

The  specific  problems  considered  are  axisymmetrically  loaded 
thick-walled  cylinders  with:  (1)  one  end  closed  with  a  flat  cap; 
(2)  one  end  closed  with  a  hemispherical  cap;  and  (3)  both  ends 
closed  with  either  flat  or  hemispherical  end  caps. 

Axisymmetric  problems  are  mathematically  two  dimensional, 
involving  the  area  which  when  rotated  about  an  axis  of  revolu¬ 
tion  generates  the  solid  of  revolution  comprising  the  region  of 
the  problem.  Since  the  geometry  of  capped  thick-walled  cylinders 
precludes  an  exact  solution,  an  approximate  numerical  solution 
was  considered. 


Of  the  numerical  methods  considered  the  finite  difference 
method  in  conjunction  with  Southwell  stress  functions  was  used1  . 
This  method  is  well  suited  to  machine  computation  and  can  be 
organized  such  that  most  of  the  tedious  labor  is  delegated  to 
the  computer.  The  axisymmetric  problem  is  essentially  solved 
when  Southwell  stress  functions  cp  and  if  are  obtained  over  the 
generating  area  of  the  solid  of  revolution  of  interest. 


Assuming  the  z  axis  to  be  an  axis  of  revolution,  the  South 
well  stress  functions  cp  and  f  over  the  generating  area  are 
governed  in  cylindrical  coordinates  by  the  equations: 


(3-la) 


ail  A  a*  _  aaco 

arar”  r  dr  +  *k.a  ”  dza 


(3- lb) 
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Stresses  for  the  axisymmetric  problem  are  defined  in  terms  of  cp 
and  f  as: 


(3-2a) 


i[S*?  [•*“-»'*] 


(3-2b) 


r 

z  r  dr 

.  lai 

zr  r  dz 


(3-2c) 


(3-2d) 


Allen  has  shown  that  boundary  conditions  are  satisfied  if  cp 
and  f  satisfy  equations  equivalent  to  the  following  on  the 
boundary  of  the  generating  area*  [See  Ref.  1,  pg.  139] 

*  -r a  (3-3a) 


to  +  (to  '  r  [t+U-v)<p]}coia  •  ro 


( 3— 3b) 


In  summary  of  the  Southwell  Stress  Functions;  Functions  cp 
and  f  must  be  found  such  that  Eqs.  (3-1)  are  satisfied  over  the 
generating  area,  and  Eqs.  (3-3)  are  satisfied  on  the  boundary  of 
the  generating  area.  Stresses  may  then  be  obtained  with  Eqs.  (3-2)  . 


134  - 


CHAPTER  2 


THE  FINITE  DIFFERENCE  METHOD  AND  SOUTHWELL  STRESS  FUNCTIONS 


Difference  Equations  in  Cylindrical  Coordinates 

Finite  Difference  Methods  are  among  the  more  frequently  en¬ 
countered  approximate  methods  for  solving  differential  equations. 
The  literature  on  this  subject  is  extensive.  Frequent  references 
are  made  in  this  study  to  the  book  by  Forsythe  and  Wasow®  . 

In  the  finite  difference  method  the  region  of  interest  is 
covered  by  a  mesh  or  grid.  The  intersections  of  the  grid  lines 
are  referred  to  as  nodal  points  and  are  classified  as  boundary 
points,  adjacent  points,  or  interior  points  according  to  proximity 
to  the  boundary.  All  differential  equations  are  approximated  by 
algebraic  difference  equations  which  involve  the  representative 
values  of  the  unknown  functions  cp  and  f  at  the  nodal  points. 
Sufficient  difference  equations  must  be  written  to  solve  for  all 
cp  and  f  values  associated  with  the  nodal  points  over  the  region. 
This  involves  approximating  Eqs.(3-1)  and  (3-3)  by  difference 
equations,  and  applying  the  equations  to  the  appropriate  nodal 
points.  Stresses  may  then  be  obtained  from  the  solved  values  of 
cp  and  f  by  difference  equations  approximating  Eqs. (3-2) . 

Where  possible,  well  known  central  difference  equations  were 
used  to  approximate  differential  equations  [See  Ref.  2,  pg.  187], 
Figure  20  shows  the  notation  employed.  Thus,  for  example, 


Figure  20.  Notation  for  nodal  points. 
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dr 


For  a  square  mesh,  which  was  used  in  the  rectangular  generat 
ing  areas,  the  difference  equations  equivalent  to  Eqs. (3-1)  are 
respectively : 

+  (1  ‘  2?)CpN  +  '"W  +  (1  +  2F>C(,S  "  =  0  <3-4a) 


t„  +  fi¬ 


ll-)* 

2r;  VN 


+  *w  +  (1  + 


2r^  ”  4*P 


CP£  +  <PW  -  2cpp 
(3-4b) 


On  the  boundary  central  difference  equations  equivalent  to 
Eqs.(3-2)  involve  nodal  points  outside  of  the  region,  i.e. 
"fictitious"  nodal  points.  Accordingly  a  backward  difference 
equation  using  the  first  three  terms  of  Newton's  backward  inter 
polating  formula  is  employed  there?  .  On  an  outer  boundary  for 
example: 

I r*k[3(  >p-4(  >s+( 

and  on  an  inner  boundary, 

h~-  k [3<  >p-4(  >N  + 


>8s]  <3-5a) 

(  )„]  (3- 5b) 


Similar  equations  can  be  written  in  the  z  direction. 


Using  Eqs.  (3-5)  and  central  difference  equations  where 
possible  to  represent  Eqs. (3-2)  yields  the  difference  equations 
for  calculating  normal  stresses  on  the  boundary  surfaces  of  a 
cylindrical  region. 


On  an  outer  surface  parallel  to  the  z  axis: 


27h  {[3  -  7l(1-v)),P  -  +  '•’aa  +  [3 


P 


-4*s  + 


(3-6a) 
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( 3— 6b) 


a  Cl  *r~  ♦  —  (l-v)lcp  -  4cp  +  cp  +  —  *  } 

0  2rh  LL  rv  J  P  S  SS  rv  PJ 

°z~  2rt  [3*P  -  4*s  +  *SS]  (3-6c) 

A  similar  equation  can  be  written  for  an  inner  boundary,  using 
the  second  of  Eqs. (3-5)  for  the  derivative  in  the  r  direction. 

On  a  boundary  parallel  to  the  r  axis: 

°r“  2*  On  -  VS  +  ♦»  -  *s]  -  7T  [*P  +  a-v)tpp]  (3-7a) 

“e25  K  -  vs]  +  P-  [*p  +  <1-''>'ppj  <3-7b> 

-  2rt[*N  -  *s]  (3-7c) 


Boundary  Conditions  for  Cylindrical  Regions 

In  general  is  is  always  possible  to  integrate  Eq. ( 3— 3a) 
around  the  boundary  of  the  region  to  obtain  if  on  the  boundary. 

To  apply  Eq. (3-3b)  it  is  necessary  to  combine  its  difference 
equation  equivalent  with  difference  equations  approximating 
Eqs.  (3-1)  in  such  a  way  that  "fictitious"  nodal  points  are 
eliminated.  Allen  discusses  in  detail  how  this  is  accomplished 
for  a  special  case  of  loading1  .  More  general  boundary  equations 
can  be  obtained  in  the  same  way.  These  equations  are  referred  to 
as  governing  cp  equations  on  the  boundary,  since  qp  must  be 
determined  there.  All  boundaries  are  assumed  parallel  to  either 
the  z  or  r  axis. 

On  an  outer  boundary, 


\ 
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On  an  inner  boundary, 


*■  - 1!1  +  r(1  +  b  (1-v)>p +  -  [3  +  7(1  +  bh 

-  rh  B»-Sf  +  U  +  77>°pr]  (3-8b> 

On  a  right  boundary, 

a  -  |jr)<PN  ♦  (i  ♦  -  %  ♦  2»„  -  2tp 

r  buvz  i 

■  r  [h  if  '  2h  °pr]  <3-8c> 

On  a  left  boundary, 

(1  -  ^)opN  ♦  (1  ♦  ^)cps  -  *p  ♦  2*b  -  2»p 

*  E  [hS  ^  +  2h  °pr]  t3-8d> 


For  cylindrical  regions  not  touching  the  z  axis  and  with 
boundaries  parallel  to  the  r  and  z  axis,  Eqs. (3-4) ,  (3-3a) ,  and 
(3-8)  are  sufficient  to  obtain  cp  and  i|r  everywhere,  providing 
the  boundary  conditions  are  in  terms  of  axi symmetric  stresses. 
However,  regions  that  have  end-caps  and  hence  touch  the  z  axis 
are  of  particular  interest  in  this  report. 


Conditions  on  the  Axis  of  Revolution 


For  solids  of  revolution  that  touch  the  axis  of  revolution 
the  governing  Eqs. (3-4)  become  undefined  in  their  present  form. 
Allen  has  pointed  out  it  is  necessary  to  apply  L'Hospital's  rule 
to  obtain  for  these  equations  on  the  z  axis  the  limiting  forms: 


dza  dza 


0 


(3-9) 
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which  implies  cp  and  f  are  linear  functions  along  the  z  axis. 
Furthermore  the  displacement  u  in  the  r  direction  relates  cp  and 
f  by  the  general  equation  [See  Ref.  1,  pg.  196] 

♦  +  (l-v)cp  =  (ff-)u 


Thus  along  the  z  axis  where  r  and  u  are  zero. 


cp 


♦ 

<l-v) 


(3-10) 


It  follows  if  i|r  is  known  on  the  left  and  right  boundaries  of  a 
region  at  the  axis  of  revolution,  then  f  and  cp  can  be  obtained 
all  along  the  axis  of  revolution  (z  axis) . 

In  the  calculation  of  t  on  the  boundary  by  integration  of 
Eq.  (3-3a)  ,  the  initial  value  of  f  at  the  starting  point  of  inte¬ 
gration  is  arbitrary.  Assume  that  the  starting  point  of  the 
integration  is  always  the  left  intersection  of  the  z-axis  with 
the  generating  area,  and  the  initial  value  of  f  there,  hereafter 

denoted  as  *  ,  is  set  to  zero.  Then  if,  denoted  by  *  ,  at  the 
L  R 

right  intersection  of  the  z  axis  with  the  generating  area  is, 

♦R  -  $-ropz  ds  s  0  <3-U) 
L 

where  the  integration  is  around  the  boundary  of  the  generating 
area  from  the  left  to  the  right  intersection  of  the  z  axis  with 
the  solid  of  revolution.  That  the  integral  vanishes  identically 
follows  from  static  equilibrium  in  the  z  direction  of  any  sector 
of  an  axi symmetrically  loaded  solid  of  revolution,  i.e., 

^apz(rA0)ds  =  0 
L 


Note  A6  is  arbitrary  because  of  axial  symmetry. 

It  may  be  concluded  from  Eqs. (3-9) ,  (3-10) ,  and  (3-11) 

that  for  *  set  equal  to  zero,  cp  and  f  will  be  zero  all  along 
L 

the  axis  of  revolution. 

With  the  values  of  cp  and  f  established  along  the  z  axis, 

Eqs. (3-4) ,  ( 3— 3a) ,  and  (3-8)  are  sufficient  to  handle  cylindrical 
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regions  that  touch  the  z  axis,  but  which  do  not  have  interior 
cavities.  Problems  that  can  be  solved  include  cylinders  with  one 
end  open,  since  they  may  be  classified  as  having  a  singly-connected 
region.  Cylinders  with  both  ends  capped  involve  doubly-connected 
regions,  and  will  be  discussed  in  Chapter  4. 


By  Eqs. (3-2a,b)  the  definitions  of  a  and  o  in  terms  of 

r  8 

and  also  become  undefined  as  r  goes  to  zero.  The  limiting 
values  of  these  equations  are: 


cp 


o 


r 


l  raa  ♦ 

a8  =  2  Ldr3 


(1+v) 


a^n 

3raJ 


Employing  symmetry  about  the  z  axis  and  letting  cp  =  iff  =  0 
on  the  z  axis  the  equivalent  difference  equation  is, 


[ 


a  1 
eJr=0 


(1+v)  (sp„>  + 


(3-12) 


Difference  Equations  in  Spherical  Coordinates 

For  cylinders  with  hemispherical  end-caps  it  is  possible  to 
employ  a  square  grid  in  the  rectangular  generating  area  and  a 
polar  coordinate  grid  in  the  quarter  ring  which  generates  the 
end-cap.  The  two  coordinate  systems  must  be  coupled  together  at 
the  juncture  of  the  end-cap  and  cylinder  as  will  be  discussed 
later. 

In  the  quarter  ring  generating  area  Eqs. (3-1) ,  (3-2),  and 
(3-3)  must  be  transformed  to  spherical  coordinates  (which  reduce 
to  polar  coordinates  for  axial  symmetry)  and  equivalent  difference 
equations  written.  Again,  central  difference  equations  were 
employed  where  possible4  .  Appendix  A  includes  the  somewhat 
lengthy  derivation  of  the  difference  equations  applying  in  the 
quarter  ring  generating  area.  The  following  equations  perform 
the  same  function  in  the  generating  rectangular  and  quarter  ring 
areas  respectively:  Eqs. (3-4)  and  (A-3) ,  (3-6)  and  (A-15)  ,  and 

Eqs. (3-8)  and  (A-12,  A-13)  . 

The  procedure  is  identical  for  the  cylindrical  and  spherical 
regions.  Of  particular  interest  is  a  problem  involving  both 
regions  coupled  together  at  the  juncture  of  the  cylinder  and  cap. 
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Juncture  of  Cylindrical  and  Hemispherical  Regions 


Figure  21  shows  the  generating  area  in  the  neighborhood  of 
the  juncture  of  a  cylinder  and  hemispherical  end  cap.  As  stated 
in  Appendix  A,  the  stress  functions  cp  and  <|i  will  be  invariant 
under  coordinate  transformation  from  physical  considerations. 
Coupling  of  the  coordinate  systems  may  thus  be  accomplished  by 
considering  the  valn.es  of  the  functions  along  the  interface  at 
points  Pi  through  r  to  be  common  to  both  regions. 


r 


In  Fig.  21,  functions  at  points  Pi  through  P4  are  assumed 
to  be  governed  by  difference  equations  in  rectangular  coordinates. 
The  horizontal  grid  lines  are  extended  into  the  polar  regions 
intersecting  radius  a-a  at  points  P2E  and  P3E.  It  is  not  neces¬ 
sary  to  extend  the  grid  lines  from  PI  or  P4,  since  the  rectangular 
difference  equations  there  involve  only  Pi,  P2  and  P3,  P4  respec¬ 
tively. 

Governing  difference  equations  equivalent  to  Eqs. (3-1)  can 
be  written  for  points  P2  and  P3  by  using  equivalent  difference 
expressions  for  the  second  derivative  in  the  z  direction  with  un¬ 
equal  intervals  [See  Ref.  1,  pg.67].  Considering  point  P3  of 
Fig.  21  for  example: 
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a3(  ) 

Jza 


]p3*  ha(l+?)  [?( 


>P3E+t 


1+ 


P3W  § 


>P3] 

(3-13) 


Using  the  above  equation  in  Eq.  (3-la)  yields  for  the  governing 
cp  equation  at  P3: 

[  (?)  (1+?t]CPP3E+  (1“lr)cpP4+  (U?)CW  (1  +  2r)cpP2“2(TS')TP3=  ° 

(3-14) 

The  value  of  the  function  at  point  P3E  in  the  above  equation 
can  be  expressed  in  terms  of  values  at  points  PP4 ,  PP3 ,  and  PP2 
by  interpolation  along  the  radius  a-a.  Using  the  first  three 
terras  of  Newton's  forward  interpolating  formula  T See  Ref.  3,  pg.  64], 

(  )p3E2!|(2-3u+ua)  (  )ow+(2u)  (2-u)  (  )  (u)  (1-u) 


PP4 


PP3 


(  *PP2] 


(3-15) 


Substituting  Eq. (3-15)  into  (3-14)  yields  the  governing  cp 
equation  at  point  P3. 

(S)  (U  [(2-3u+ua)cppp4+(2u)  (2-u)cppp3-(u)  d-u)cppp2] 

+  (1-  27)8PP4  +  -  21 ’  0 


(3-16) 

A  similar  equation  can  be  written  at  point  P2  by  using  the 
appropriate  nodal  points  and  u  and  §  values. 

The  governing  equations  for  f  at  point  P3  can  be  found  by 
combining  the  operators  of  Eqs.  (3-1  )  and  (3-13)  . 

76)  (iVD  [2-w  >  W (2u)  (2-u)  *PP3  -  (»)(!-«)  *ppj 

+  (1"  2r)*P4+<l+i')*P3w+<1+  2r)*P2  "  2  (-5^)  *p3 

[<2-3u+u’)cppp4+(2u)  (2-u)cppp3  -  u(l-u)cppp2] 


?d+?) 
2 


^i+;^p3w  +  V^ps 


=  o 


(3-17) 
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A  similar  governing  f  equation  can  be  written  for  point  P2, 
again  using  the  appropriate  nodal  points  and  values  of  u  and  £. 

The  functions  along  the  radius  a-a  at  points  PPl,  PP2 ,  etc. 
are  assumed  to  be  governed  by  appropriate  equations  developed  in 
Appendix  A  for  spherical  regions.  The  difference  equations  at 
those  points  will  involve  points  Pi  through  P4  without  any  modi¬ 
fications  being  necessary.  Likewise  points  P2W  and  P3W  will  be 
governed  by  unmodified  rectangular  difference  equations  and  will 
involve  points  P2  and  P3.  Thus  a  sufficient  number  of  equations 
are  available  to  solve  for  all  unknown  values  of  cp  and  h  in  the 
neighborhood  of  the  juncture  of  the  cylinder  and  hemispherical 
cap. 


143  - 


CHAPTER  3 


APPLICATIONS  TO  CYLINDERS  WITH  ONE  END  CAPPED 


Cylinder  with  a  Flat  Cap 

Figure  22  shows  the  generating  area  for  an  example  problem 
with  a  flat  cap.  The  loading  considered  was  a  unit  normal  pres¬ 
sure  p  on  the  outer  surface.  Normal  stresses  obtained  on  the 
boundaries  by  the  finite  difference  method  and  Southwell  stress 
functions  are  also  shown. 

An  indication  of  the  error  in  the  solution  on  the  boundaries 
parallel  to  the  z  axis  can  be  obtained  by  comparing  the  or  values 

calculated  from  the  cp  and  if  functions  to  the  known  applied  normal 
boundary  stress.  The  error  is  obviously  greatest  at  the  re¬ 
entrant  corner,  where  the  boundary  was  rounded  off  to  avoid  a 
singularity.  The  nodal  point  there  is  considered  to  be  an  interior 
node,  but  the  error  inherent  in  representing  a  rapidly  changing 
geometry  by  a  set  of  discrete  points  is  apparent.  Grading  the 
mesh  around  the  corner  would  hopefully  reduce  the  error  there, 
but  would  require  more  hand  labor  and  machine  storage.  A  stress 
concentration  obviously  exists  at  the  corner  which  cannot  be 
found  directly  by  the  finite  difference  method,  but  hopefully 
could  be  approached  with  a  closer  spacing  of  nodal  points.  The 
solution  approaches  the  Lame'  solution  quite  rapidly  as  z  in¬ 
creases  beyond  the  neighborhood  of  the  end  cap. 

In  obtaining  the  solution  for  the  example  problem  computer 
programs  were  designed  and  used  for:  (1)  generating  most  of  the 
coefficients  in  the  matrix  of  coefficients  for  the  unknown  cp  and 
t|i  values  at  interior  nodal  points;  (2)  solving  by  an  iteration 
procedure  for  cp  and  (r  from  the  complete  matrix  of  coefficients 
and  constants  vector;  and  (3)  calculating  normal  stresses  on 
boundaries  parallel  to  the  z  axis.  These  programs  are  included 
in  Appendix  B  Program  B-l  employs  Eqs.(3-4),  while  program  B-3 
employs  Eqs. (3-6) . 

vhe  generation  of  the  coefficients  for  boundary  points  of  a 
rectangular  generating  area  must  be  done  by  han.*  using  Eqs.  (3-8)  . 

It  is  also  necessary  to  modify  the  coefficients  in  the  equations 
for  the  points  adjacent  to  the  boundary  by  inserting  appropriate 
boundary  coefficients  or  constants  according  to  Eqs. (3-4)  .  The 
constants  vector  is  derived  from  applying  Eqs. (3-8)  and  (3-4). 
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Figure  22.  Stresses  in  an  open  cylinder  with  flat  cap. 
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Figure  23.  Hoop  stress  on  the  outer  surface  of  a  cylinder 


The  iteration  procedure  for  solving  the  difference  equations  and 
the  program  for  calculating  normal  stresses  on  the  boundary  may 
be  used  without  hand  modification. 

In  the  second  example  problem,  a  unit  ring  of  pressure  was 
applied  at  a  particular  location.  Figure  23  shows  the  hoop 
stress  on  the  outside  boundary.  The  main  objective  of  this  ex¬ 
ample  was  to  study  how  rapidly  the  stress  decays  away  from  the 
point  of  load  application.  The  stress  apparently  decays  quite 
rapidly,  implying  the  effect  of  the  end  caps  also  decays  rapidly 
as  z  increases  beyond  the  neighborhood  of  the  end. 


Cylinders  with  a  Hemispherical  Cap 

Figure  24  shows  the  generating  area  for  an  example  problem 
with  hemispherical  cap.  The  loading  considered  was  a  unit 
pressure  on  the  outer  cylindrical  portion  of  the  solid  of  revolu¬ 
tion.  Normal  stresses  on  the  boundaries  are  also  shown. 

Again  an  indication  of  the  error  involved  may  be  obtained 
by  comparing  the  calculated  value  of  stress  normal  to  the  boundary 
to  the  applied  stress.  As  might  be  anticipated,  the  maximum  error 
occurs  near  the  juncture  of  the  two  regions.  Again,  as  in  the 
case  of  a  re-entrant  corner,  grading  the  mesh  near  the  juncture 
would  probably  reduce  the  error  there.  Additional  study  of  the 
coupling  of  two  coordinate  systems  for  a  finite  difference 
solution  would  be  desirable. 

In  obtaining  the  solution  for  the  cylinder  with  a  hemi¬ 
spherical  end  cap  the  three  computer  programs  described  in  the 
previous  sub-section  were  used.  In  addition  programs  were 
written  and  used  for:  (4)  generating  the  coefficients  in  the 
matrix  of  coefficients  for  all  points  in  the  quarter-ring 
generating  area;  and  (5)  calculating  normal  stresses  on  the 
xnner  and  outer  surfaces  of  the  hemispherical  cap.  These  pro¬ 
grams  are  included  in  Appendix  B.  Program  B-4  employs  Eqs. (A-5)  , 
(A-6) ,  (A— 12)  and  (A-13) .  Program  B-5  employs  Eqs. (A- 15)  and 
(A- 16)  . 

Generation  of  the  matrix  of  coefficients  corresponding  to 
Fig.  24  requires  the  use  of  programs  B-l  and  B-4.  In  addition, 
for  nodal  points  on  and  adjacent  to  the  boundaries  of  the  rectan¬ 
gular  generating  area  the  procedure  described  in  the  last  sub¬ 
article  which  involves  hand  modification  must  be  used.  For  the 
interior  nodal  points  along  the  interface  of  the  two  regions 
equations  similar  to  Eqs. (3-16)  and  (3-17)  must  be  applied  by 


147 


-2p 


148 


Figure  24.  Stresses  in  a  cylinder  with  hemispherical  cap. 
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hand.  The  final  matrix  of  coefficients  must  be  in  a  form  accept¬ 
able  to  the  iteration  procedure  program,  B-2.  Obviously  there 
should  be  no  duplication  of  numbering  of  unknowns  over  the  two 
regions  comprising  the  total  generating  area.  Programs  B-3  and 
B-5  will  scan  the  total  output  of  the  solution  by  program  B-2, 
and  calculate  the  appropriate  stresses. 

The  hand  labor  involved  for  this  problem  comes  about  mainly 
from  hand  applications  of  difference  equations  to  the  boundary 
points  of  the  rectangular  generating  area  and  to  points  along 
the  interface  between  the  two  areas.  The  time  expended  in  machine 
calculations  is  spent  mainly  in  obtaining  the  solutions  to  the 
difference  equations. 
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CHAPTER  4 


ANALYSIS  OF  CYLINDERS  WITH  BOTH  ENDS  CAPPED 


Cavities  in  Axisymmetr ^.c  Solids 


Allen  has  discussed  the  solution  for  ax i symmetric  solids  with 
a  cavity  [see  Ref.  1,  pp.  195-198].  He  shows  that  for  any  axisym- 
metric  body  single  valuedness  of  displacement,  v,  in  the  z  direc¬ 
tion  is  insured  by  the  equation 


+ 


+  l^jjdr  dz  *  0 


0-18) 


where  C  is  any  closed  curve  in  the  generating  area.  For  any 
curve  not  enclosing  a  cavity,  the  above  integral  equation  re¬ 
quirement  is  equivalent  to  Eqs.(3-1).  However  for  solids  with 
cavities  not  touching  the  axio  of  revolution,  the  above  equation 
applied  around  the  inner  boundary  can  be  used  to  determine  the 
no  longer  arbitrary  initial  value  of  f  there.  For  if  the  initial 
value  of  f  on  che  exterior  boundary  is  selected  to  be  arbitrary, 
then  the  initial  value  of  i|i  on  the  inner  boundary  is  not. 

For  an  axisymmetric  solid  with  a  cavity  touching  the  axis  of 
revolution  the  above  equation  yields  no  information  about  ♦  when 
applied  around  the  inner  boundary.  In  this  case  it  is  not  possible 
to  roturn  to  the  starting  point  of  integration  because  of  the  "gap" 
along  the  z  axis.  If  an  integration  is  performed  in  the  negative 
r  region,  then  the  integral  all  around  the  axisymmetric  cavity 
correctly  vanishes  identically  since  this  is  essentially  a  re¬ 
tracing  in  the  oposite  direction  of  the  integral  in  the  generating 
area  above  the  z  axis.  Another  method  must  be  used  to  determine 
the  initial  value  of  ♦  on  an  inner  cavity. 


Determination  of  t  on  an  Inner  Boundary 

Consider  the  solid  of  revolution  shown  in  Fig.  26.  From  a 
previous  discussion  in  chapter  2  the  values  of  f  and  cp  will  be 
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related  constants  along  the  z  axis.  The  behavior  of  f  on  an  inner 
boundary,  denoted  by  f  along  r  *  p,  as  p  approaches  zero,  is  of 

particular  interest. 
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Figure  26.  Determination  of  i|i  on  the  inner 
surface  of  a  closed  cylinder. 

It  will  be  argued  that  <i  approaches  the  value  of  i|i  on  the  z 

axis,  denoted  by  ^  ,  as  p  approaches  zero.  Thus,  the  initial  value 

z 

of  f  on  an  inner  boundary  at  the  axis  of  revolution  is  equal  to 
the  arbitrary  initial  value  of  if  on  an  outer  boundary  at  the  axis 
of  revolution. 

Consider  the  physical  effect  of  making  p  ever  smaller.  As  p 

decreases  o  in  the  central  "core"  of  material  along  the  z  axis 
z 

will  increase,  but  at  a  decreasing  rate.  For  as  p  is  made  smaller 
the  central  core  will  carry  less  of  the  load  that  ties  the  ends 
together.  This  force  is  distributed  to  the  cylinder  walls,  which 
are  theoretically  capable  of  carrying  any  load.  Furthermore  the 
elongation  of  the  core  in  the  z  direction  will  remain  small,  and 
will  become  essentially  constant  as  p  approaches  zero.  The  latter 
statement  is  true  because  the  displacement  of  the  cylinder  ends  in 
the  direction  of  the  z  axis  is  affected  to  an  ever  smaller  extent 
by  the  decreasing  of  p  after  the  cylinder  walls  are  carrying  almost 
all  of  the  load  formerly  carried  by  the  core  of  the  cylinder.  It 
follows  the  strain  and  hence  the  stress  on  the  core  for  small  values 
of  p  remains  defined  and  essentially  constant  as  p  becomes  smaller. 
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Recall  that 


c 

z 


I  £1 

r  dr' 


and  along  the  z  axis, 

da 


°z  =  ■  3r»  • 


Thus  if  c?z  remains  defined  throughout  the  core  as  p  becomes 
small,  then  -j-  throughout  the  core  and  —Jf-  on  the  z  axis  remain 

defined.  This  implies  that  f  approaches  t  as  p  approaches  zero. 

X  z 

It  is  therefore  possible  to  set  f  =  0  on  the  inner  boundary 
at  the  z  axis,  if  the  initial  value  of  f  on  the  external  boundary 
at  the  z  axis  is  also  set  to  zero.  If  the  cylinder  is  subject 
to  internal  loads  on  the  inner  cavity,  then  f  will  vary  there 
according  to  Eq. (3-3a) . 

The  initial  value  of  Jr  on  the  interior  surface  of  a  thick 
walled  cylinder  with  both  ends  capped  may  be  found  by  the  above 
analysis.  The  solution  otherwise  is  identical  to  the  solution 
for  a  cylinder  with  an  open  end. 

Figure  27  shows  normal  stresses  in  a  closed  flat-end  cylinder. 
The  cylinder  was  loaded  by  a  uniform  pressure  over  the  ends.  The 
analysis  discussed  above  was  used  in  the  determination  of  the 
initial  value  of  f  on  the  inner  boundary  of  the  cylinder.  Values 
of  o  and  a  were  not  plotted  at  the  reentrant  corner  as  it  was 

IT  Z 

felt  the  values  obtained  there  were  unrealistic  due  to  the  rapidly 
changing  geometry  and  the  relatively  coarse  mesh  employed. 
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CHAPTER  5 


i 

SUMMARY 

A  method  has  been  presented  for  the  determination  of  stresses 
in  axially  loaded  thick  walled  cylinders  with  one  end  capped.  The 
end  cap  may  be  either  flat  or  hemispherical.  Computer  programs 
have  been  designed  to  handle  most  of  the  tedious  labor  associated 
with  the  finite  difference  technique.  Example  problems  of  cylinders 
with  both  types  of  end  caps  are  presented. 

The  maximum  error  in  the  example  problems  considered  occurs 
in  the  neighborhood  of  reentrant  corners  and  at  the  juncture  of 
the  cylinders  and  end  caps.  Additional  refinement  of  the  grid 
would  probably  reduce  the  error  in  these  locations,  but  would 
require  additional  hand  work  to  obtain  the  corresponding  difference 
equations. 

An  analysis  of  thick  walled  cylinders  with  both  ends  capped 
has  been  made.  It  is  shown  that  on  the  axis  of  revolution  the 
initial  value  of  the  Southwell  stress  function  f  on  an  inner  sur¬ 
face  may  be  set  equal  to  f  on  the  outer  surface.  An  example 
problem  employing  the  above  analysis  is  presented  for  a  cylinder 
with  both  ends  capped  by  flat  caps. 

Of  additional  interest  would  be?  (1)  the  effect  of  re¬ 
fining  the  grid  at  the  reentrant  corners.  (2)  a  more  effective 
coupling  together  of  the  two  coordinate  systems  at  the  juncture 
of  the  cylinders  and  hemispherical  end  caps,  and  (3)  cylinders 
with  end  caps  that  are  not  complete  hemispheres. 


i 

I 

( 
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APPENDIX  A 


DEVELOPMENT  OP  FINITE  DIFFERENCE  EQUATIONS  FOR 
SOUTHWELL  STRESS  FUNCTIONS  IN  SPHERICAL  REGIONS 


As  noted  previously  axi symmetric  elasticity  problems  are 
mathematically  two  dimensional  and  involve  the  generating  area 
o£  a  solid  of  revolution.  Thus  for  spherical  regions  the  prob¬ 
lem  involves  a  generating  area  which  can  most  readily  be  described 
by  the  coordinates  (R,@)  as  shown  in  Fig.  28. 


z 


Figure  28.  Coordinate  system  and  notation. 


Since  the  differential  equations  governing  the  Southwell 
stress  functions  are  of  second  order  it  is  feasible  to  transform 
them  into  the  coordinates  (R,0)  and  use  appropriate  polar  coor¬ 
dinate  finite  difference  representations.  From  physical  con¬ 
siderations  the  stress  functions  f  and  cp  will  be  invariant  under 
coordinate  transformations.  The  equations  defining  the  stresses 

o  ,  a c  and  t  become  in  the  (R,0)  coordinates: 
r  0  z  zr 


c 


r 


1 _ 

R  Sin8 


d  ($+♦ )  Cos0 
8R  +  R 


te±tl  _ 1 _ 

d©  R  Sin© 


[♦+  (l-v)rpjj 


c 


e 


V _ 

R  Sin® 


{[ 


Sin® 


dtp  Cos®  dcp~| 

&r  r  aej 


l _ 

vR  Sin© 


(a) 


(A-l) 


[♦  (l+v)cpjj 


(b) 
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1 _ 

R  Sin@ 


(A-l) 


i - r 

R  Sin@  L 


cos©  ajL  _  sine  ii-i 
S R  R  S®J 


The  normal  stresses  a  and  oQ  can  be  found  by  stress  trans- 

R  W 

formations  of  the  type  involving  direction  Cosines  1,  m,  and  n: 


a  *  a  la  +  a  ma  +  21m  t 
R  z  r  zr 


Substituting  Eqs. (A-l)  in  the  stress  transformation  equations 

yields  a  and  art. 

R  © 

aD  38  p ' c‘:" %  {sin©  (“!“)  -  Cp,S~^  (—^“)  +  Sin3©  ("[Sin© 

R  R  Sin©  l  SR  R  S©  L  SR 


(A- 2) 


©  R 


Cos3© 

R  Sin© 


[♦+  (l-v)cpj}. 


The  governing  equations  for  the  Southwell  Stress  Functions 
cp  and  i|r  in  the  coordinates  (R,®)  are: 


Sacp  1 


[cos®  ft  -  =  °- 


ill  _  1_  fCot8  it  .  ilil  _ii£eCosa@.2ifi2_ 

9Ra  r3  LCot0  30  je!J  IT08  0 


3cp  Sin©  Cos© 

MR - 5 -  (A_3) 


dcp  Sin3 0  rdcp  Sin©  Cos  ©  Sacp  Sin3 0 
+  SR  R  +  d©  Ra  +  S0a  Ra  * 

(b) 
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For  a  circular  boundary,  ,  a  *  (90° -0).  The  boundary 

O  M  O  R 

conditions  for  the  Southwell  Stress  Functions  may  then  be  written 
in  terms  of  (R,0)  as: 


and 


(a) 


{||sin®  + 


5cp  Cos0 
50  R 


RSIM  [♦+d-v)ep]}sin@ 


(A-4) 


=  R  Sin0  a 

pr 


(b) 


The  resulting  problem  in  a  hollow  hemispherical  region  can  be 
solved  by  finite  difference  techniques  in  a  fashion  similar  to 
that  for  the  hollow  cylindrical  region.  Functions  cp  and  f  must 
be  found  that  satisfy  Eqs. (A- 3)  over  the  region  and  the  boundary 
conditions  (A-4)  on  the  boundary  of  the  region.  To  accomplish 
this  central  difference  equations  in  polar  coordinates  will  be 
employed.4 

The  notation  shown  in  Fig.  29  will  be  used  in  the  develop¬ 
ment  of  finite  difference  equations  in  the  coordinates  (R,0)  . 


Figure  29.  Notation  for  nodal  points  in  polar  coordinates. 
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Only  equal  intervale  in  the  radial  direction  will  be  considered. 
The  incremental  angles  11  shall  likewise  be  equal.  The  "compass" 
notation  is  used  to  identify  the  unknown  values  of  cp  and  f  at  the 
nodal  points  labeled  as  shown  above. 

The  finite  difference  equations  corresponding  to  Eqs. (A-3a) 
and  (A-3b)  are  respectively: 


and 


,n  cot©-)  ,n  coten  ^r.  ,n  .  *1 

+  V  >N +  ww +  V  L1^— >S-2[1+W  >p=° 

(A- 5) 


*E+ 

-  [Cosi®  +  h  sin* ®]cP E  +  [lh?in®  COS0>„E-  (^)3sin® 


[sin0+TK:ose]cpN  -  [sin®  Cos®  Jp^-fcos^-jjSin3®]^ 
(2riR'[Sin®  Cos®>sW  “  W  Sin©  [sin®  “  nCos©]cps 


(2^R)[Sin®  Cos®]cP se  =  0 


(A-6) 


As  might  be  anticipated  the  boundary  conditions  are  more  com 
plex  in  the  (R,@)  coordinates  than  in  the  (r,z)  coordinates.  The 
equation  (A-4a)  may  be  integrated  around  the  boundary  to  obtain  i|t 
on  the  boundary,  as  was  done  for  the  (r,z)  coordinates.  The  equa 
tion  (A-4b)  has  for  its  finite  difference  equivalent: 

li  h 

Sina@  cp  +  (— )  Sin©  Cos©  cp  -  Sina®  cp  -  (—r)Sin©  Cos©  rp( 
E  Y^R  N  W  T]R  I 

-  (-jp)  (1-V)cpp  +  tE-i |!w-  =  2hRSin®  (A-7) 

The  "east"  or  "west"  values  will  be  regarded  as  "fictitious" 
values  to  be  eliminated  according  to  whether  the  equation  is 
applied  on  the  outer  or  inner  boundary  of  the  region. 


Consider  for  example  the  devel upment  of  the  boundary  equation 
governing  cp  on  the  outer  boundary  of  the  polar  generating  area. 
Note  again  that  i|(  is  known  on  the  boundaries  from  Eq.  (A-4a)  . 
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To  apply  Eq.  (A-7)  on  the  outer  boundary  the  cp  and  *  terms 

E  E 

must  be  eliminated  by  combining  Eqs.(A-5),  (A— 6)  and  (A-7). 

Equation  (A-5)  is  multiplied  by  (-Sin9®)  and  added  to  Eq.  (A-7) 
to  yield  an  equation  free  of  cp_,  i.e. 

£-Sina0  (—)  (l-r|  -C-j— )  +  (“)  Sin©  Cos9]cpN-2^SinaeJcpw 
^-Sina®(~)  (1-t-ri  C^t— )  -  (~)Sin®  Cos®Jcps+2^Sina® 

[1+(^),]-  R(1-V)>P  +  -  »W  -  r  *P  *  2hR  Sin0  °pr 

(A-8) 


To  eliminate  i|!  from  the  above  Eq.  (A-8)  it  is  necessary  to 

E 

obtain  a  second  equation  free  of  cp  .  A  combination  of  Eqs.  (A-5) 

E 

and  a  modified  form  of  Eq. (A-6)  will  yield  an  appropriate  equation 

The  modified  form  of  Eq. (A-6)  employs  a  backward  difference 

representation  of  the  cross  derivative  of  (p  in  Eq. (A- 3b)  .  This 

is  necessary  to  eliminate  the  fictitious  values  of  cpVT_  and  cp 

NE  SE 

in  Eq.  (A-6) ,  for  there  is  no  possibility  of  otherwise  handling 
these  terms.  The  backward  difference  formulas  used  on  the  outer 
and  inner  boundaries  for  the  cross  derivative  were  respectively: 


dacp 

d0dR 


?.atP 

dOdR 


2t|h  [^N  "  ^NW  +  ^SW  “  ^sl 


2rih  [^NE  “  +  ^S  ^Se] 


(A- 9) 


Using  the  first  of  Eqs. (A-9)  in  the  representation  of  Eq. 
(A-3b)  yields: 

*E  +  <V[i-^>N+v(V[1+^]v2[1+4>a>P 


-[c°sa0  +  2iSina®}pE  -  (^)Sin®  [cos0(--  1)+ 


-  (“j)  [sin®  CosQjcp^  -[cosa9-~Sina0jcpw+(^)[sin©  Cos®]cp 


sw 
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+  (^)  Sine[co.e  (|  -  1)  -  ^Sine}ps  -  0 


(A-10) 


Multiplying  Eq.  (A-5)  by  £cosa0  +  —SinaPj  and  adding  to  the 

above  equation  yields  a  second  equation  free  of  cp  . 

£ 

+  'pE+{(^)a[1-2Sfii][CO8a0+25Sin,0]'[^Sin@][Co8e(|-  11 


+  ^isin0]K  -  <Msin@  CO80>iro+(I)[sina0},w  +  (J^> 

[sine  cose]cpsw  +  {(^>’[1  ♦  ■QSr£][cosae  + 

+  Kr51"®!0080  (t  ‘  11  ‘  ^R®in0]K  =  0  (A'U) 

Subtracting  Eq. (A-8)  from  Eq.  (A-ll)  yields  a  boundary  equation 
for  cp  on  an  outer  boundary  free  of  all  fictitious  values. 

(“~)  {[ l-^p]!; Her811*3 ®]~[ Sina e]"[ TlSinQ  Cos@J|  c 

-  (^>[Sin0  C089>NW+  (Sir'a0)[2+R'}(,W+  (^R>[Sin®  CO80>SW 
+  ^ [  It j [ 1-K^gS ina  © ]~ [~ S ina  9 ]•>■  [ “nS in©  Cos@JJ<ps 

*2fe>a{1-  [<f>a  <t>  a-v)]-Sin>8[l-(^)a4«f)a+l]]K 


Note  that  in  the  above  equations  f  f  ,  and  )  will  be  known 

NS  P 

values  on  the  boundary.  The  above  equation  is  the  difference 
equation  for  cpp  on  an  outer  boundary,  R  *  constant. 
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In  a  similar  fashion  the  governing  equation  for  cpp  on  an  inner 

boundary  may  be  obtained.  It  is  necessary  in  this  case  to  use  the 
second  of  Eqs.  (A-9)  in  the  difference  representation  of  the  cross 
derivative  in  Eq. (A- 3b) .  Thi  resulting  governing  cp  equation  is: 

[sin’0(2-£)}pE  ♦  <^) [Sine  (Jj)' ‘ {[l  -  ***■] 

[l-^Sin»eJ-[sinae]-[TiSine  CoaeJ}«N+  (^)  [[l  + 


[l-^Sina9]-[sin*®J+[tiSine  Co»9 J jtpg-  [sine  Cos0]cpgE 

-2(^>S{1+t  <1-n-sina9[i-(f)a4«f)a  * 


*  2hR  Sine  o 

pr 

In  the  above  equation  (as  in  A-12)  the  values  of  *  ,  and 

N  S 

*  are  know.. 

TP 

If  only  hollow  hemispherical  regions  joined  to  cylindrical 
regions  are  studied  as  in  this  report,  then  Eqs. (A-5)  ,  (A-6) , 
(A-12)  ,  and  (A-13)  are  sufficient  to  obtain  cp  and  f  over  the  cir¬ 
cular  area  generating  the  hemispherical  region.  The  value  of  cp 
will  be  known  along  0  =  0  in  general,  and  the  values  of  cp  along 
0  *  tt/2  are  assumed  to  be  determined  by  appropriate  difference 
equations  written  for  the  adjoining  rectangular  area  generating 
the  cylindrical  region.  More  detail  on  this  final  statement  is 
given  in  the  main  body  of  the  report. 

The  values  of  cp  and  f  obtained  by  the  above  equations  may 

then  be  substituted  into  the  difference  equation  representation 

of  Eqs. (A-lo,d)  and  Eqs. (A- 2)  to  yield  values  of  o0 ,  t  ,  o 

9  zr 

and  o  . 

Since  the  values  of  cp  and  f  are  solved  on  the  boundary  and 
in  the  interior  of  the  region,  it  is  necessary  to  use  a  backward 

difference  formula  for  .  The  first  three  terms  of  Newton's 

dR 

backward  interpolating  formula  for  the  first  derivative  may  be 
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used  here.3  On  the  outer  boundary  the  radial  derivative  operator 
is  represented  by: 


aj _ i 

3R 


k[  3<  >p  - 4(  >w+  < 


and  on  the  inner  boundary, 


3R  2h  L 


3(  )p  -  4( 


>E  + 


(A- 14) 


The  double  subscript  implies  the  value  of  the  function  at 
the  next  node  beyond  the  west  or  east  nodal  point  in  the  radial 
direction. 


Using  Eqs.  (A- 14)  and  central  difference  equations  for  the 
angular  derivatives  in  Eqs. (A- lb)  and  (A-2)  yields  the  difference 
equation  relations  for  normal  stress  on  an  outer  boundary. 


V 

2hR 


{ (i*)  [coteKw  (^>  [cote>s+[3^iwi> 


[vRSina®]*p} 


(a) 


0R*  *2hR  {Sina®[(^)[CotS]'pHW,MW-%-(?i)[COt0]t|)S 

+[  3' S]»p]- 4)[COt8]*H+W4*W  <A‘15) 

+  (^)rcot0]»s  +  [3 (b) 

aaet  2hR  {C08S®[(^)[COte],,,N^WM-4<pW-<^)tCOte^S 

+  [  3‘rii^e]cpp]+  V  [cote]v  w4*w 

-  (^)[cote]»g  -  [3  +  j*<cote)»]:p}  (c) 

On  an  inner  boundary  the  difference  equations  which  are  used 
to  determine  the  normal  stresses  are: 
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°6*  *2hR  {HpEE+4V  (^>  [COteK-  (^»  [COt9>s 


r-i  2(1"v)  i  . 

r  2h  1 

L3  vRSina0  J^P  + 

LvRSina0j 

°R9£  2hS  {Sin30[^EE+4q)E^  (^tCOt8l'PN-(^i)CCOte]'PS 

(A-16) 

+  r^0]cpp]^EE+4*E“(^)[COt0]*N 
+  (^)Ccot8]*s.[3  +  f>p}  (b) 

"0*  {CO88e[^EE+4cpE+(^R)[COt0]cpN  ‘  (^)[C°t0]*S 

+  RSina0^  top]  +^EE"4^E+(riR^Cot0-'^N 

-(^Ccotells  +  t3  -  |^(cote)a]tp}  (c) 
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FOREWORD  TO  APPENDIX  B 
Computer  Programs 


The  programs  described  were  written  in  Fortran  II  for  the 
IBM  1620.  In  certain  programs,  such  as  the  coefficients  generator 
program  B-l,  the  use  of  cards  is  implied.  It  is  convenient  to 
remove  certain  machine  generated  "flagged"  cards  by  sorting  and  to 
replace  them  with  hand  written  cards  according  to  the  special 
situation  encountered. 

Descriptions  of  programs,  flow  charts,  and  source  deck  list¬ 
ings  are  included.  The  stress  calculation  program,  B-5,  applying 
uO  hemispherical  end  caps  is  highly  specialized  and  serves  primarily 
as  an  example  program  which  can  be  modified  for  problems  other  than 
the  example  solution  shown  in  Fig.  24.  A  table  of  contents  for 
Appendix  B  follows  for  convenience  in  locating  a  specific  program. 


Program  Page 

B-l.  Coefficients  Generator  for  Southwell  Stress 

Functions  in  Rectangular  Coordinates  173 

B-2 .  Equations-solver  Program  by  Iteration  185 

B-3 .  Calculation  of  Normal  Stresses  on  Boundaries 

Parallel  to  the  Axis  of  Revolution  197 

B-4.  Coefficients  Generator  for  Southwell  Stress 

Functions  in  Polar  Coordinates  209 

B-5.  Calculation  of  Normal  Stresses  on  Boundaries 

of  a  Hemispherical  Cap  224 
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APPENDIX  B 


COMPUTER  PROGRAMS 


Program  1. — Coefficients  Generator  for  Southwell  Stress  Functions 
in  Rectangular  (r,z)  Coordinates. 

This  program  essentially  generates  a  matrix  of  coefficients 
corresponding  to  linear  algebraic  equations  derived  by  finite 
difference  techniques  and  the  use  of  Southwell  Stress  Functions. 

It  generates  the  matrix  of  coefficients  of  the  unknown  values  of 
cp  and  f  occuring  at  interior  nodal  points  of  a  square  grid  cover¬ 
ing  the  regions  of  interest.  For  a  complete  solution  of  any  prob¬ 
lem  the  coefficients  generated  here  must  be  supplemented  by  coef¬ 
ficients  obtained  by  hand  for  boundary  points  and  interior  points 
adjacent  to  the  boundary.  This  program  will  flag  unknown  boundary 
values  affecting  nodal  points  adjacent  to  the  boundary  as  .JWl, 

JW2,  JEl ,  etc.,  which  must  then  be  inserted  by  hand.  The  result¬ 
ing  output  is  in  a  form  suitable  for  input  to  the  iteration  program 
included  in  this  study. 


Definitions  of  input  parameters  which  describe  the  region  of 
interest  and  identify  the  unknowns  follow. 


K  *  last  row  number  in  the  region,  obtained  as  distance 
measured  in  grid  intervals  h  from  the  z  axis  to  the 
furtherest  parallel  grid  line  in  the  region. 


J  **  number  associated  with  the  initial  unknown  cp  value 
which  is  assumed  to  occur  at  the  left  end  of  the 
initial  row.  All  unknowns  are  numbered  sequentially 
from  left  to  right  and  in  the  increasing  r  direction, 
through  the  cp  and  then  if  unknowns. 


rinit 


first  row  number  in  the  region,  obtained  as  distance 
between  the  z  axis  and  the  nearest  parallel  grid  line. 


Limits  of  Rows  *  distances  measured  in  grid  intervals  h 
from  a  common  r  axis  to  the  first  and  last  nodal 
points  respectively  of  rows  first  through  last. 


Input—  (fixed  point  throuohout) 

1.  Parameter  Card 

Cols  1-5  K 

Cols  6-10  J 

Cols  11-15  I. 

mit 
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2.  Limits  of  Rows  Cards — (Rows  I  ^  through  K  in  sequence,  one 

card  per  row) 

Cols  1-5  Row  No. 

Cols  6-10  Distance  to  first  nodal  point  in  row 

Cols  11-15  Distance  to  last  nodal  point  in  row 


Output 


Cols  1-10 
Cols  11-15 
Cols  16-20 


Coefficients  of  Unknowns 
Rows  of  Unknowns 
Cols  of  Unknowns 


An  additional  comment  should  be  made  that  does  not  deal 
directly  with  this  program.  Before  writing  by  hand  the  equa¬ 
tions  involving  the  boundary  nodal  points  the  map  of  the  grid 
covering  the  region  should  be  drawn  to  avoid  error.  Note  the 
unknowns  on  the  boundaries  parallel  to  the  z  axis  should  be 
numbered  from  left  to  right,  as  this  is  required  by  the  stress 
calculation  program  B-3. 
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COEFFICIENTS  GENERATOR  IN  RECTANGULAR  COORDINATES 
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SET  SWITCH  FOR  2ND  VARIABLE  AND  INITIALIZE 
IS  VI 

KA-IRFf  UNIT) 

Jl*Jl2*Hl 
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07  FORMAT! 1 15. 2X, 3HJNI 
00  FORMAT! II5.2X.3HJW2 
09  FORMAT! I 15»?XV 3HJE2 
END 


Program  2 .--Equations-Solver  Program  by  Iteration 

The  program  described  in  the  following  pages  involves  an 
iteration  procedure  which  reduces  to  Gauss  Seidel  iteration  if 
the  incorporated  over-relaxation  factor  is  set  to  one.  This  was 
done  for  most  of  the  problems  of  this  study,  as  the  equations 
resulting  from  the  finite  difference  method  in  conjunction  with 
Southwell  stress  functions  tended  to  diverge  or  show  little 
benefit  from  the  use  of  an  over-relaxation  factor  other  than 
one. 


The  program  is  described  in  terms  of  general  x  unknowns. 

For  this  study  the  x  unknowns  represent  the  Southwell  stress 
functions  cp  and  i|r.  Thus  two  unknowns  exist  at  each  interior 
nodal  point  of  a  region  and  one  at  each  boundary  nodal  point. 

Immediately  following  the  source  deck  listing  of  this  pro¬ 
gram  is  a  brief  conversion  program  that  converts  the  output  from 
any  cycle  of  iteration  to  a  form  suitable  for  input  of  initial 
"guesses"  of  the  unknowns  for  a  following  cycle  of  iteration. 

In  this  way  the  solution  can  be  interrupted  at  the  end  of  any 
cycle  of  iteration  and  then  continued  at  a  later  time. 
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Young-Frankel  Iteration  Program  for  Solution  of 
Sparsely  Populated  Matrices. 


The  purpose  of  this  program  is  to  solve  a  very  large  number 
of  simultaneous  equations  by  the  Young-Frankel  over-relaxation 
procedure  [See  ref.  2,  p.  242].  The  size  of  the  system  that  may 
be  processed  is  determined  only  by  the  number  of  non-zero  coef¬ 
ficients  and  constants.  In  solving  various  types  of  differential 
equations  by  finite  difference  methods,  a  system  of  linear  equa¬ 
tions  is  obtained.  The  resulting  matrix  of  coefficient  terms  is 
usually  sparsely  populated  and  the  distribution  may  or  may  not 
follow  a  pattern.  Ordinary  elimination  techniques  are  impractical 
where  the  coefficients  are  randomly  distributed,  because  a  large 
amount  of  intermediate  computer  storage  may  be  required.  Itera¬ 
tive  techniques  are  usually  more  time-consuming,  but  the  storage 
requirements  may  be  kept  to  a  minimum.  Instead  of  storing  the 
coefficient  array  in  its  usual  two-dimensional  array,  where  most 
of  the  elements  are  zero,  two  one-dimensional  arrays  are  stored; 
one  is  the  list  of  non-zero  coefficients  and  constants  arranged 
in  order,  the  second  is  an  index  matrix  giving  the  relative  posi¬ 
tion  of  the  corresponding  element  in  the  coefficient  and  constant 
matrix.  The  ordering  of  the  elements  of  both  matrices  can  be 
illustrated  in  the  example  shown  below: 


Order 

>11*1+  >la*s 

®aa*l+®aa*a*® 

*51*1  +-»»»X»’Cs 

List  of  Coefficients 

Index  Matrix 

1 

an 

1 

2 

al  3 

3 

3 

*88 

5 

4 

*83 

6 

5 

*11 

7 

6 

*33 

9 

7 

Ci 

1 

8 

c3 

3 

For  a  system  of  n  equations,  an  additional  set  of  w  loca¬ 
tions  is  used  for  intermediate  approximations  of  the  unknowns, 
which  can  be  set  at  the  discretion  of  the  operator.  As  each 
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new  x  is  computed,  the  next  approximation  is  taken  according 
to  the  following  equation,  which  forms  the  basis  of  the  Young- 
Frankel  method: 

x,(l,+1)  *  */*)  +  r,  (x/-x,  (*)) 

where 

x^k+l )  a  next  approximation  to  xt  (corrected) 
x/k )  =  preceeding  approximation 

th 

r,  =  relaxation  factor  for  the  i  equation 

x,'  =  next  approximation  to  xt  (uncorrected) 

Note  that  if  the  relaxation  factor  =  1,  then 

x(k+*)  =  xt' 

which  is  the  Gauss-Seidel  procedure.  It  has  been  demonstrated 
that  the  rate  of  convergence  may  be  greatly  accelerated  by 
choosing  a  factor  other  than  1  in  the  correction  step.  A  defi¬ 
nite  optimum  may  exist  which  will  give  the  maximum  convergence 
rate.  For  a  given  problem,  however,  this  factor  must  be  deter¬ 
mined,  usually  by  trial  based  on  previous  experience.  The  pro¬ 
gram  is  arranged  to  allow  a  different  relaxation  factor  for 
each  equation,  if  desired. 

Additional  features  of  the  program  include  multiple¬ 
processing  of  constant  vectors  without  re-entry  of  the  coef¬ 
ficient  list,  and  intermediate  output  of  results  which  may  be 
obtained  at  any  time  during  processing.  Since  initial  approxima¬ 
tions  may  be  entered,  the  program  may  be  stopped  at  any  time  by 
output  of  intermediate  results  which  may  then  be  re-entered  as 
approximations.  A  check  routine  is  included  to  ascertain  whether 
all  diagonal  elements  of  the  coefficients  matrix  have  been 
entered.  A  maximum  value  of  the  total  allowable  relative  error 
in  the  solution  vector  may  be  specified;  at  any  time  that  the  sum 
is  less  than  that  given  the  program  will  give  the  results  and 
halt.  A  maximum  number  of  iterations  may  also  be  given;  when 
this  number  has  been  completed,  the  program  stops  and  the  results 
to  that  point  are  given  as  output. 


Input 


Type: 

1.  Parameter  card 
Cols  1-10 


TOLR,  maximum  allowable  relative  error 
(fixed  point) 
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Cols  11-15  No.  of  equations  (fixed  point) 

Cols  16-20  Maximum  No.  of  iterations  (fixed  point) 

Cols  21-48  Problem  identification  (any  alphameric 

data) 

2.  Coefficient  cards  (one  card  per  element) 

Cols  1-10  Element  (floating  point) 

Cols  11-15  Row  No.  (fixed  point) 

Cols  16-20  Column  No.  (fixed  point) 

(No  zero  coefficient  is  permitted  as  input) 

(The  row  and  column  numbers  are  understood  to  give  the 
position  of  the  element  in  the  original  two-dimensional  array) 

(Cards  must  be  sorted  in  order  Cols  11-20) 

3.  One  blank  card. 

4.  Relaxation  factors  (as  many  cards  as  desired— minimum  of 
one,  maximum  of  n) 

Cols  1-10  Factor 

Cols  11-15  Highest  row  No.  for  which  this  factor 

is  to  be  used. 

If  one  factor  is  used  then  the  Cols  11-15  must  contain 
No.  of  equation. 

If  more  than  one  factor  is  used,  then  these  cards  must 
be  sorted  in  order  by  Cols  11-15. 

5.  One  blank  card. 

6.  Initial  approximations  (as  many  cards  as  desired — minimum 
none,  maximum  of  n — all  others  assumed  zero) 

Cols  1-10  Initial  approximation 

Cols  11-15  Position  of  corresponding  x  element  in 

solution  vector. 

Cards  must  be  sorted  by  Cols  11-15. 

7.  One  blank  card. 
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Young-Frankel  Iteration  Program 
for  Solution  of  Sparsely  Populated  Matrices 
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Increase  maximum 
allowable  itera¬ 
tions  by  present 
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140  CALL  SETUP  00098 

SOLVE  BY  ITERATION  00099 

160  CALL  SOLVE  00100 

FINAL  PART  00101 

TEST  RESIDUAL  SUM  AND  NO.  OF  ITERATIONS  00102 


IF(RSUM-TOLR) 170*100* 180  0010) 

170  CALL  OUT  00104 

GO  TO  240  00108 

180  CALL  SSWTCHf 1* J) 


{ 


o  **  eg  e*  «  ^  o  **  «#  tr  seo>o  •*  m«*  eg  m  *  ir\  <o  * 

******  ******  oj  eg  nnn  eg  eg  eg  m  m  moooooon 

******  **********  ******  *4  *****  «*  -*  **000000** 
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ooo  ooooo  ooo  oooo  o  oooooooo 
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IS*- I  00135 

D0330K«1*N  00136 

IF 1 1 V-I A) 300* 310v300  00137 

300  IFfMI J)-K)310»3209310  00138 

320  AIISN)-A(J)  00139 


*  4  4  4  *  *  *  V>OOOOOOlAIAIMMMfMMAlA«  4  0  0  0  0  0 
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540  IEME+1  00201 

550  IFI IE-I I ) 560» 560* 570  00202 

560  J*JM  00203 

G0T0500  00204 

GET  NEM  X  AND  OELTA  X  00205 


O  O  O  O  •«  *4  **  •* 

NNNNNNNNN  N  N  N  N  N  N  N  M 

ooooooooo  oooooeoo 
ooooooooo  oooooooo 


Converter  Program 


1  READ  100 ,I,A,J,B,K,C,L,D 

100  FORMAT (4 (14 , E13 . 5 , 3X) ) 
PUNCH101,  A,I,B,  J,C,K,D,L 

101  FORMAT (F10. 5, 15) 

GO  TO  1 

END 
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Program  3. — Calculation  of  Normal  Stresaes  on  Boundaries  Parallel 
to  the  Axis  of  Revolution. 

This  program  calculates  normal  stresses  on  boundaries  paral¬ 
lel  to  the  z  axis  according  to  Eqs.  (3-6) .  It  employs  the  output 
from  any  cycle  of  iteration  of  the  equation  solver  Program  B-2. 

To  use  the  output  of  Program  B-2  the  identifying  preliminary  cards 
must  be  removed  from  the  particular  output  being  used. 

The  equations  (3-6)  are  applied  from  left  to  right  along  the 
boundaries.  This  program  assumes  unknowns  will  be  numbered  se¬ 
quentially  from  left  to  right,  however  there  may  be  gaps  in  num¬ 
bering  between  rows.  The  Equation  Generator,  B-l,  numbers  in¬ 
terior  points  from  left  to  right,  but  care  must  be  taken  to  num¬ 
ber  boundary  points  also  from  left  to  right.  Known  values  of  Jr 
from  left  to  right  on  the  boundaries  must  be  read  in  sequentially 
after  each  data  card.  If  Jr  is  constant  on  the  boundary,  then 
only  one  Jr  card  is  required  and  the  sense  switch  1  must  be  turned 
on. 


Definitions  of  Input  Terms. 

KEQUA  =  number  of  equations  in  input  data  from  program  B-2. 

PR  =  Poisson's  ratio  of  material 

KEEP  =  identifying  number  of  first  unknown  at  left  end  of 
boundary.  See  Fig.  30. 

L2  =  identifying  number  of  unknown  cp  value  at  first  interior 
point.  See  Fig.  30. 

L3  =  identifying  number  of  unknown  Jr  value  at  first  interior 
point.  See  Fig.  30. 

L4  =  identifying  number  of  unknown  cp  value  at  second  interior 
nodal  point. 

L5  =  identifying  number  of  unknown  Jr  value  at  second  interior 
nodal  point. 

LAST  *  identifying  number  of  last  unknown  at  right  end  of 
boundary. 

R  *  distance  to  boundary  in  terms  of  grid  spacing  H. 

H  *  grid  spacing.  This  number  is  1.0  for  all  work  done  in 

this  study.  In  general  it  depends  on  the  units  in  which 
the  original  difference  equations  are  written. 
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Figure  30.  Notation  for  program  3. 


SIGN  *  +1.0  for  outer  boundary  calculations,  -1.0  for  inner 
boundary  calculations. 

PSI  *  value  of  the  second  Southwell  stress  function  t  at  boundary 
points  at  which  normal  stresses  are  calculated. 

PSIC  *  constant  value  of  f  at  boundary  points. 


Input 

1 .  Header  Card 

Cols  1-4  KEQUA,  Right  Justified  Integer  (Fixed  Pt.) 

Cols  11-20  PR  ,  Floating  Point 

2.  Solutions  Cards — Output  from  Program  B-2. 


3. 


4. 


Data  Cards 


Cols:  1-3 

KEEP 

4-6 

L2 

7-9 

L3 

10-12 

L4 

13-15 

L5 

16-18 

LAST 

19-28 

R 

29-38 

H 

39-48 

SIGN 

PSIC  Card  (S.S. 

1  on)  , 

,  Right  Justified  Integer 
,  Right  Justified  Integer 
,  Right  Justified  Integer 
,  Right  Justified  Integer 
,  Right  Justified  Integer 
,  Right  Justified  Integer 
,  Floating  Point 
,  Floating  Point 
,  Floating  Point 

or  PSI  Cards  (S.S.  1  off) 


Cols  1-10  PSIC  or  PSI 


(Fixed 

(Fixed 

(Fixed 

(Fixed 

(Fixed 

(Fixed 


Pt.) 

Pt.) 

Pt.) 

Pt.) 

Pt.) 

Pt.) 
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Output 

Sigma  R 
Sigma  Theta 
Sigma  z 


values  of  normal 
values  of  normal 
values  of  normal 


stress  at  point 
stress  at  point 
stress  at  point 


in  r  direction, 
in  6  direction, 
in  z  direction. 
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Calculation  of  Normal  Stresses  on  Boundaries 
of  a  Thick-walled  Cylinder 


INPUT 
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L5M.541 

PUNCH  10»(IvSIGR(I)»SIGT<I)«SIGZtI)« IsKEEP«LAST) 

FORMAT! I4,3X,3E16.8,25X) 

FORMAT! 52HP01 NT  SIGMA  R  SIGMA  THETA  SIGMA  Z) 


Program  4.— Coefficients  Generator  for  Southwell  Stress  Functions 
in  Polar  (R,6)  Coordinates. 


This  program  generates  the  matrix  of  coefficients  correspond¬ 
ing  to  the  linear  algebraic  equations  derived  by  finite  difference 
techniques  and  Southwell  stress  functions.  It  essentially  applies 
Eqs.(A-13),  (A-5) ,  (A-12) ,  and  (A-6)  to  a  quarter-ring  which,  when 
rotated,  generates  a  hemispherical  solid. 

With  reference  to  Figure  31,  this  program  will  generate  the 
entire  matrix  of  coefficients  for  the  generating  ring  area,  except 
for  the  coefficients  of  the  unknowns  along  6  -  n/2.  These  unknowns 
lie  on  the  interface  between  the  end  cap  and  cylinder,  and  equa¬ 
tions  applying  here  must  be  written  by  hand  according  to  Eqs.(3-16) 
and  (3-17) .  Note  however  that  this  program  will  write  the  complete 
equations  for  points  adjacent  to  the  interface,  and  therefore  the 
unknowns  along  the  interface  must  be  numbered.  There  must  be  no 
duplication  in  numbering  of  unknowns  anywhere  over  the  ring  and 
subsequently  coupled  rectangular  area.  The  output  of  this  program 
is  in  a  form  similar  to  that  of  Program  B-l,  and  can  be  combined 
with  it  as  input  to  the  equation  solver  Program  B-2. 

Referring  to  Figure  31  the  equations  of  Appendix  A  are  applied 
in  the  following  orders  (1)  Eq.  (A-13)  on  the  inner  boundary  start¬ 
ing  at  point  "MIKE"  and  moving  in  the  increasing  0  direction, 

(2)  Eq.  (A-5)  to  interior  points  in  the  increasing  0  and  then  R 


r 


Figure  31.  Notation  for  program  4. 
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directions,  (3)  Bq. (A-12)  on  the  outer  boundary  in  the  increasing 
8  direction,  and  (4)  Eq. (A-6)  to  interior  points  in  the  increasing 
0  and  then  R  directions.  The  first  three  steps  involve  ep  unknowns 
at  the  inner  boundary,  interior,  and  outer  boundary  nodal  points 
respectively.  The  last  step  involves  f  unknowns  at  interior  nodal 
points. 

Definitions  of  Input  parameters  and  symbols  used  in  the  pro¬ 
gram  follow. 

IA  ■  Interior  radius,  an  integer  in  terms  of  h,  the  grid 

spacing  in  the  R  direction,  >  0 

IB  *  Exterior  radium,  an  integer  in  terms  of  h,  >  IA+1. 

IJK  ■  Number  of  angular  subdivisions  of  one  quadrant,  23. 

MIKE  ■  Identifying  number  associated  with  the  unknown  cp 

value  on  the  inner  boundary  at  the  initial  point  of 
the  ring.  See  Figure  31. 

MAC  -  Identifying  number  associated  with  the  unknown  cp 

value  on  the  inner  boundary  at  the  initial  point  of 
the  juncture  of  the  ring  with  the  cylinder.  See 
Figure  3. 

PR  *  (1-v) ,  one  minus  Poisson's  ratio  of  the  material. 

PHI  »  Input  values  of  cp  along  the  inner  and  outer  boundaries. 

SIGPR  ■  o  *  Input  values  of  applied  stress  in  r  direction 

pr  along  inner  and  outer  boundaries. 

ETA  *  Increment  of  angle  in  polar  grid. 

THETA  ■  Polar  angle  at  location  of  equation  application. 

See  Fig.  31. 

Xl,X2,etc.  ■  Varying  coefficients  of  cp  and  if  values  appearing  in 
Eqs. (A-13) ,  (A- 5) ,  etc. 

CONST  *  Constant  coefficients  of  cp  and  if  values  appearing  in 

Eqs. (A-13) ,  (A- 5) ,  etc. 

CONH  *  Constant  term  of  an  equation  that  appears  in  the  con¬ 

stants  vector.  It  is  made  up  of  input  boundary  values 
of  f  and  boundary  stress  o  according  to  Eqs. (A-13), 
(A-5) ,  etc.  pr 


Input 

1.  First  Data  Card 
Cols  1-3 
Cols  4-6 
Cols  7-9 


IA  A  right  justified  integer  (Fixed  Pt.) 
IB  A  right  justified  integer  (Fixed  Pt.) 
IJK  A  right  justified  integer  (Fixed  Pt.) 
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Cola  10-12  MIKE  A  right  justified  integer  (Fixed  Pt.) 

Cola  13-15  MAC  A  right  juatified  integer  (Fixed  Pt.) 

Cola  16-30  PR  Floating  Point. 

2.  Additional  Data  Carda — Use  floating  point  numbers,  8  per  card, 

10  cols,  per  number,  blank  carda  read  aa  zeroa.  Data  ia  sequen¬ 
tial,  and  no  spaces  must  be  left  unless  they  represent  zeroa. 

Example  Input  data  is  listed: 

Inner  Boundary  Terms: 

^[mike]'  a[pr  mike]'  ^[mike+1]'  a[pr  mike+1] '  '  ^[mike+IJK-2] , 

Outer  Boundary  Terms:  (follow  above  data  with  no  spaces) 

*[mike+  (IJK-1)  (IB-IA)  ] '  q[pr  (same  pt.]'  *[mike+  (IJK-1)  (IB-IA)  +  1] 

a[pr  same  pt.]'  '  ^[mike+  (IJK-1)  (IB-IA+l)-!]  '  g[pr  siune] 

Horizontal  Boundary  Terms:  (follow  above  data  with  no  spaces) 

^[mike]  '  ^[mike]'  ^[mike+IJL] '  ^[mike+IJL] '  ^[mike+2  (IJL)  ]  ' 

same] ' ’ * ’ '  ^mike* (IJK-1) (IB-IA)],  *[same] ' 

Lower  Point  on  Vertical  Boundary:  (follows  above  data  with  no 
spaces) 

^[mike+ IJK-1] 

Top  Point  on  Vertical  Boundary:  (follows  above  data  with  no  spaces) 
^[mike+  (IB-IA+1)  (IJK-1)] 


Output 

Coefficients  Cards 

Cols  1-10  Coefficients  in  Floating  Point 
Cols  11-15  Row  number  of  Coefficients 
Cols  16-20  Column  Number  of  Coefficients 


Constants  Cards 

Cols  1-10  Constants  in  Floating  Point 

Cols  11-15  Corresponding  Row  Number  of  Constant 
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Coefficients  Generator  in  Polar  Coordinates 


1 


Calculate 
MN,  CONH,  Ml 

* 

Punch  XI,  I,  MN) 
X4,  Ml,  CONH 


I>K 


I<K 


<5 


Calculate 
Ml,  CONST 


•© 


II-L 


0—3 


Initialize 

IQ-0 

MZ-IA+II 

PII-M2 

HRHO-l./PH 

n-kk+ijn 


Calculate 

MN-MAC+1 

I*KK 

IQ-IQfl 


w 

f  Punch  X3,  X6,\ 

^3l\— 

k  1/  MAC,  MN 

< 

© 


Punch  X6,  I,  MN, 
X2,  X5,  Ml,  CONST 


I<KJ 


I>KJ 


I-I+l 


17 


Initialize 

NN*»2 

KK-MIKE+IJL 
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Calculate 
X10,X11 ,X12 ,Ml 
CONST 


OUTPUT 


[Punch  CONST, 
ll,  Ml,  Xll 


I<KK 

I>KK 


Calculate 
CONR,  MN 

r4 - 

' Punch 
I,  XI 

CONH,' 

2.  MN  1 

I>N 


KK«N+l| 


Punch  X12 , 
I ,  MN 


I<N 


I-I+l 


IKIE 


II-II+l 


-  20  k  - 
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EQUATION  GENERATOR, FOR  ONE  QUADRANT  OF  A  CYLINDRICAL  SECTION 
IJK'THE  NO.  OF  DIVISIONS  OF  ONE  QUADRANT, MUST  BE  AT  LEAST  3. 
THERE  MUST  BE  AT  LEAST  ONE  INTERIOR  LINE. 
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Program  5. — Calculation  of  Normal  Stresses  on  Boundaries  of  a 
Hemispherical  End  Cap. 

This  program  calculates  values  of  normal  stress  on  the  boun¬ 
daries  of  hemispherical  regions  according  to  Eqs. (A-15)  and  (A-16) . 
As  the  program  is  now  written  it  is  specialized  for  the  generating 
quarter  ring  area  shown  in  Fig.  24.  Thus  r\  =  tt/30  in  the  equa¬ 
tions  (A-15)  and  (A-16) .  However  this  program  could  be  easily 
modified  to  a  more  general  form. 

Definitions  of  input  parameters  follow: 

L  Number  identifying  first  unknown  in  the  output  of 

program  B-2  (lowest  row  No.  in  the  matrix  of  coeffi¬ 
cients  for  cp  and  )  . 

N  Number  identifying  last  unknown  in  output  of  program 

B— 2 . 

K  Integer  multiplying  “n  such  that  0=K“n,  see  Fig.  32. 

JJ=2  Defines  location  in  terms  of  K  of  initial  point  at 

which  stress  is  calculated  on  inner  boundary. 

IK=13  Defines  location  in  terms  of  K  of  final  point  at 

which  stress  is  calculated  on  inner  boundary. 

MM  Subscript  of  cp  at  IK  from  map  of  region.  (See  Fig. 

32)  . 

MK  Subscript  of  cp  at  initial  point  at  which  stress  is 

calculated  on  the  outer  boundary. 

NK  Subscript  of  cp  at  final  point  at  which  stress  is 

calculated  on  the  outer  boundary 

II=L+70=71  Subscript  of  cp  at  initial  point  on  inner  boundary 

for  the  example  problem  shovvn  in  Fig.  24.  This  num¬ 
ber  should  be  changed  for  a  different  numbering  of 
the  initial  point.  This  requires  a  modification  of 
the  program,  i.e.,  number  of  70  must  be  changed. 

KK=N-78=82  Subscript  of  cp  at  point  preceding  final  point  on 
inner  boundary  at  which  stress  is  calculated. 

Special  case  for  Fig.  3-5,  number  of  78  must  be 
changed  for  different  problem. 
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Input — Fixed  Point 


1.  First  Card 


Cols.  1-5 

L 

Cols.  6-10 

N 

Cols.  11-15 

K 

Cols.  16-20 

JJ 

Cols.  21-25 

IX 

Cols.  26-30 

MM 

Cols.  31-35 

MX 

Cols.  36-40 

NX 

2 .  Data  Cards 

Output  of  Program  B-2,  all  preliminary  cards  removed. 

3 .  Output 

Values  of  Sigma  Theta  *  o  ,  Sigma  Nor  *  c_»  Sigma  Tan  -  n  , 
versus  K.  Note  Kt)  *  0.  9 


Figure  32.  Notation  for  program  5. 
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Calculation  of  Normal  Stresses 
on  Boundaries  of  a  Hemispherical  End-Cap 
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